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The Use of Exponential Sums in Step by Step 
Integration 


Introduction. In the step by step numerical open integration of a system 
of ordinary differential equations 


(1) fila, Sn, 

filzi, «++, Zn, #) is regarded as a function of t, F(t), whose values are known 
tt+h 

at t,t — h,---,t — (m — 1)h and for which f F(r) dr is desired. It is 
t 


customary to assume that F(t) is a polynomial of degree » — 1 whose 
coefficients are determined by the given values of F, i.e., the values at 
t,t —h,---,t — (wm — 1)h.1 If the polynomial F(#) is specified by these 
conditions, then the desired integral is a linear combination of the values 
F(t), F(t — h), ---, F(t — (m — 1)h). 

GREENWOOD’ has pointed out that the polynomial assumption for F(?) 
may not always be the most desirable one. In particular it might be assumed 
that F(t) is a linear sum of exponentials and still obtain that the integral is 
a linear combination of the values F(t), F(t — h), ---, F(t — (m — 1)h) with 
coefficients independent of ¢. The authors also were led to this conclusion 
by their own previous work.? 

Sections 1-5 constitute a discussion of practical procedures and contain 
a description of the results obtained. Sections 6-9 contain the proofs of 
the results. The. present discussion is for an individual step. The final total 
error effect on the solution and stability considerations is not given. 

1. In an open integration procedure, the expression 


(2) h(aoF(t) + ai F(t — h) + + — (m — 1)h)) 
is used to approximate i F(r) dr. The error e(#) is given by 


(3) e(t) = A(aoF(t) + a F(t — h) + + — (m — 1)h)) 


tt+h 
= f F(r) dr. 
t 
This error function is linear in F(#) and hence if it is zero for F(t) = e’*; 
4 = 1, ---,m, it will be zero for any F that is a linear combination of these 


exponentials. The condition that e(t) be zero for F = e’* can be written in 
the form 


when one cancels out e”* and divides by h. Letting x; = e~’*, one obtains 
the equation 


do + + + + = — (1 — x)/{x In x}. 


These equations for different »; will yield m linearly independent equa- 
tions in do, @;, «++, @,-; since the determinant of this system is the Vander- 
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mondian 
1, 


Notice that the a; depend only on the product »,h. 
In the step by step integration procedure an open integration is used first 
to obtain from the function values at t,t — h, ---,t — (m — 1)h, a set of 


values of the unknowns 2; and f; at t = t + h. One uses the new values at 
th 

t+ h in a “closed” integration formula to recompute F(r) dr. While 
t 


this method may have certain disadvantages in general, the authors have 
employed alternate open and closed steps to obtain a simple estimate of 
the truncation error with a minimum of computation.* In a closed inte- 
gration ~~ F(t + h) is assumed to be known and the approximating func- 


tion to F(r) dr is h[boF(t + h) + F(t) + + — — 2)h)]. 
Thus the error of the closed procedure e,(#) may be written 
(3c) ec(t) = h[boF(t + h) + F(t) + + — (m — 2)h)] 


t+h 
_ f F(r) dr. 
t 
The 5; are chosen under the assumption that the error is zero for 
F(t) = e’#;4 = 1,2, ---, . This leads to the system of equations 
(4c) bo + + = (1 — 


Unless otherwise stated the subsequent discussion refers to the open 
integration procedure. 

2. One may wish to consider the same value »; more than once in the 
enumeration of the »;. For instance, one might want a certain frequency 7 
to appear three times among the »;, i.e., one might set v) = » = vm = v3. 
It will be indicated below that this will yield superior accuracy for fre- 
quencies near ». Higher order contacts for various frequencies may be 
obtained as follows. Suppose in equation (3), F(#) = e. Let the error e(t) 
be expressed in the form he(Ahk)e™. Then by certain obvious cancellations, 
equation (3) yields 


(6) €(Ah) = do + ayem™™ + + — (e™ — 1) (AK). 
The previous discussion established = for = v2, Now 
for a double » it is clear that e(vk) = 0 and — eon), = 0. e(Ah) may 


be considered a function of « where u = dh and “ = a e(u) = he’(u). 


Thus the partial derivative condition is equivalent to de(u)/du = 0. 
Other functions of u may be used as independent variables and here it 
is quite convenient to introduce x = e~*. Then 


(6’) €(x) = do + + + + f(x) 
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where 


(7) f(x) = (1 — x) In 


The condition de/du = 0 is equivalent to de/dx = 0, and the system of 
equations de/du = 0, d*e/du? = 0 is equivalent to de/dx = 0 and d*e/dx*? = 0. 
Similarly if m = » = ve = v3, let x; = e*. The first three equations 
e(x;) = 0, e(x2) = 0, e(xs) = 0 of type (6’) may be replaced by the three 
equations e(x) = 0, de/dx = 0, d’e/dx? = 0, at x = x,. Hence, we may use 
the three equations 


+ + + = — f(x), 
(8) Q, + + + — = — f'(x), 
+ Gage; + + — 1)(m — = — f" (x1), 


to determine the a; instead of the first three equations of (4’). It is clear 
how multiplicities in frequencies other than x; can be taken care of by 
similar methods. One can readily show that the determinant of the new 
system is not zero. 

The use of the equations (8) is generally the most convenient way to 
determine the a; when they are subject to a ew condition. The 
conditions 


de 
e(u) = 0 


=0,--- for «= 0, 


are readily seen to be precisely the condition that e(#) = 0 for F = 1, t, #, 

-,?”—. Thus the polynomial case is included in the above discussion 
provided multiplicities are considered. The final results of this paper are 
not affected by letting a number of the frequencies coincide and from the 
practical point of view, multiple frequencies may be used whenever it is 
convenient. 

3. The foregoing method of integration can be effectively applied to 
functions of the form 


(9) F(t) = + off) 


where a(t) is small but otherwise arbitrary. From equation (3) it is seen 
that the error made in integrating the function F by the above method is 


(10) e(t, F) = hE + elt, 0), 


the significance of the double argument being clear. ¢(A;k) is given by 
equation (6). 

The error may be estimated by equation (10) provided «(Ak) for 
= «++, A, and e(t, are studied. The terms of the summation are then 
readily estimated since |e“| is given by exp (¢ Re (A)). e(Ak) is evaluated 
as follows. Suppose that the integration method is precise for step 4 and 
complex frequencies ---, Now let 


y=e™M—il, —1, 
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Then 

(11) = TT G- y) |(4. + nti + + 
j= 


where 


(12) Skat = yyy 
+ 


and the A; are given by the series 
(13) y((1 + y) In (1 + y)) = Ao + Ary + + 
(Cf. Theorem I of §7 below.) The A; may readily be calculated since 


If the integral is expanded by the binomial theorem and integrated term 
by term, one obtains® 


(-1)" 


A,= r(r + 1)(r + 2) (7 + (nm — 1)) dr. 
The first eleven values of A; are listed below. 

A, = 1 Ae = .31559 
A, = —.5 Az = —.30422 
A> = .41667 Ag = .29487 
A; = —.37500 Ay = —.28698 
A, = .34861 Ay = .28019 
As = —.32986 


The A; decrease monotonically in absolute value, and alternate in sign. 
An expression for an overestimate of 


A, + + Ani2S2, n41 cee 


can readily be given. Suppose there is a Y,0 < Y < 1, such that |y| < Y 
and |y;| < Y for 7 = 1,2,---,#. Then Lemma 2 below shows that 


n 
|An + A n41 + n41 + | < |An|(1 + (n + 1)Y 

+ (n+ 1m + + +++) = — 
Thus the error is 


(14) < [i ly - yy, 


However this overestimate may be inconvenient when Y is rather close 
to unity, in which case the original formula (11) will probably give a much 
smaller estimate. 


We note the effect of the factor [J |y — y;|. Since y — yj; = e** — e~*4, 
j=1 


if 
tl 
ti 
“4 
fe 
=) 
x oO 
h 
ti 
if 
C 
e 
il 
re 
| 
Ww 
= al 
ce 
re) 
a 
5 
hi 
re 
V; 
| 
of 
cl 
Vj 
el 
m 
n 


rm 


ign. 


| 


n+1) 


lose 
uch 


—vih 
’ 


USE OF EXPONENTIAL SUMS IN STEP BY STEP INTEGRATION 67 
if Xk and v;h are small, this difference is similar to h|\ — »;| ; consequently 
this factor is essentially 4" [JT | — »;|. The original objective of this inves- 

j=l 


tigation was to establish that the error goes to zero like h**+'. Otherwise a 
favorable comparison with the known results on polynomial integration 
could not be established. 

It should be emphasized that the above results hold even when a number 
of the »; and consequently y; are equal. For a » which occurs r times, we 
have a factor |\ — vo|* in II, which indicates a small ¢ for A near y. 

On the other hand one can also take care of the case in which a par- 
ticular \4; has multiple weight in the expression (9) for F. For instance, 
if \ has double weight, one finds a term Cte“ in (9) replacing some term 
C;e*. This new term Cte“ can be written 0(Ce™)/d\ and consequently in 
equation (10) instead of Cje(A;h)e', there appears 


8(Ce(Ah)e™) /AX. 


Higher weights for \ will correspond to higher partials with respect to 
in the obvious manner. 

A similar investigation of the closed integration procedure yields the 
result 


where the S;,n4:, y, yi are defined above. B; = A; + Ais, By = Ao. (Cf. 
Corollary of Theorem I below.) 

The important properties of closed integration are present in this result. It 
can be shown by a gross approximation that form > 1, |B,| < |A,|/(m — 1) 
and B,A, <0. For m — 1 211, the error of the closed integration pro- 
cedure is materially smaller than the error of the open integration and of 
opposite sign. Two purposes are thus achieved by following an open inte- 
gration by a closed integration. The error is reduced and a good step trun- 
cation overestimate may be found. 

The error effect due to the function ¢ is calculated in Theorem II. An 
overestimate similar to those above is given in the Corollary. 

4. The use contemplated for this type of integration is the solution of 
a system (1) in which the functions f;(z, ---, Z,, 4), when obtained as func- 
tions of time, F(t) are of the form (9). In such a use, the complex frequencies 
d; are not known until the numerical work has been done. On the other 
hand, one may have initial estimates of reasonable accuracy and the above 
results indicate that one might well use these initial estimates either as the 
v; directly or to indicate a region in which the »; can be chosen to determine 
the constants a; used in the integration. This means of course that the set 
of »; must be chosen in advance of the integration. One hopes that the 
choice of the set of »; is such that for each frequency A in the F; there is a 
v; such that the factor |\ — »;| will be small. If the set of \ has many more 
elements than the number of »; permitted by practical considerations, a 
multiple »y; may have to compensate for a number of A. 

In iterative procedures, an open integration is usé¢d*to yield the first step 
values of the unknown z at ¢ + h. These z are used to compute the f at 


= 


68 USE OF EXPONENTIAL SUMS IN STEP BY STEP INTEGRATION 


t+ h and these values of f are then used to obtain a new value of the 
increment in z for the interval (t,t + 4). The resulting new values of z 
may be used to recompute the f. This process may be repeated until no 
significant change occurs. But the result even then is not perfect.’ If the f; 
are complicated, the amount of work involved may be proportional to the 
number of iteration steps and better accuracy may be available from a 
procedure of simple open integration doing the same amount of work for 
the same interval of the independent variable but using a smaller step h. 
Thus, for example, an open step followed by a closed step for h = At may 
yield a result less accurate than two open steps of length 4 = At/2. 

On the other hand, in order to use the closed integration procedure to 
estimate the truncation error, i.e., to take advantage of the fact that in 
general the error for the closed integration has the opposite sign to that of 
the open and much smaller, one need not recompute the f. The f obtained 
by substituting the results of the open integration from ¢ to ¢ + h can be 
used to obtain the Az increment for the closed step from ¢ to t + h and also 
the increment for the open integration from ¢ + h to t+ 2h. After the 
comparison between the increments Az from the open step and the closed 
step for the interval ¢ to ¢ + h, the closed increment is ignored and the open 
process is continued. 

It also should be pointed out that the true open truncation error is not 
given by equation (3) under these circumstances but by 


tt+h 


where the f is from equation (1) and %, ---,2%,, are true solutions of the 
system of differential equations which agree with the computed solutions 
at t. However, in general, for a fixed frequency, 7(#) may be obtained from 
e(t, e) by means of a factor which is close to one when u = Ah is small. 
A similar situation exists for the closed integration truncation error. (Cf. 
references in footnote 3.) The stability of these integration procedures 
should also be considered. 

5. In the appendices to this paper appears certain supplementary mate- 
rial of an illustrative nature or designed to assist the application of the 
results. In view of formula (11), above, one would normally require that 
|y| «1 for every frequency present so that ¢(Az) will be small. This 


yields an upper bound for the h’s which one can use. A quick way of finding _ 


this upper bound is indicated in Appendix 1. Since the determinant of the 
system of equations (4) is small, a straightforward numerical elimination 
will involve a great loss of accuracy in the computation of the a;. Practical 
formulas for the accurate calculation of the a; are given in Appendix 2 with 
a method of successive improvement. Appendices 3 and 4 contain numerical 
results which illustrate the theory. 

6. In order to evaluate effectively the error it will be necessary to con- 
sider briefly the higher Vandermondians, i.e., determinants of the form 


1, “1, x1", 
1, X2, x2", 


(15) Vilar, Xn) = 


1, Xn, Xn2, 


7 
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To study this we must consider the result of taking successive divided 
differences of powers. 
We define 


Jn = — x1")/(x2 — x) 
LemMA 1. For r < n, 


(16) [x, le = > eee 
---+ar=n—r+1 


(in this summation a; = 0,1,2,---). Forr =n +1, 


(16’) [x1, = 1 
and forr >n + 1, 
(16’’) = 0. 


For r = 2, we see immediately that 


In general, if the formula for r < n is true for r — 1, we have 


n—r+2 
= 2 x 
++ 
+ 


Sele 


n—r+2 
a=1 agt ++ 


Every monomial obtained by expanding this expression has 
Bi + By + bony =n—r+i, 


and conversely, every possible term of this type occurs. Consequently for 
r<n 
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To prove the remaining statements, let r = m. Then 
[x1, Sa — + Xe + eee + Xn 


and 

The result for r > m + 1 is obvious. 

The definition of [x1, ---,x,], can be extended to the cases of multiple 
x; by means of equations (16), (16’), and (16”). 

LemMA 2. For r<n-+1, the number of terms in the expansion of 
[x1, «++, is the binomial coefficient (, 

Subject to the condition, a; + --- + a, =m —r-+ 1, one has as many 
ways of choosing a, ---, a, as one can arrange » — r + 1 different things 
in r groups, blank groups being admissible. This last is the specified bino- 
mial coefficient. 

Now the higher Vandermondians may be evaluated. 


1, Xo, +, Xe™—2, 
1, 
0, [x1, x2 hh, [x1, X2 |n—14k 
= II («i — . 
0, [m, [x, Xn |n—2, [x, Xn 


This can readily be seen by subtracting the first row of the original deter- 
minant from each of the others and then factoring out x; — x, from the 
ith row for i = 2, ---, 2. One next subtracts the second row from the third, 
fourth, ---, #th row and again factors out x; — x2 for i > 2. Continuing 
this process, one will obtain the displayed result. In the second determinant 
the terms below the diagonal are zero and all but the corner term of the 
main diagonal are unity. Thus 
Lemma 3. 


The determinant obtained from a Vandermondian by differentiating a 
row and setting the variable of this row equal to the variable of another 
row will be considered. More generally, a determinant may have the same 
variable in a number of rows each of which, except the uppermost, is 
obtained by differentiating the one above it. For these a notation is intro- 
duced which may be illustrated. 


Vi( {xi}, {x2}, ee +) 


0, 1, 2m, — 2)x:***, (n — 1+ 
0,0, 2, +++, (m —2)(m — (n —1+ k)(m — 24+ 
0, 1, 2x2, (m — 2)x2"-*, (n — 1+ 
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It is clear, for instance, that 


Vil {xr}, Xa, Xn) Vi(x1, x.) | 
and that 


Vo(x1, =0 
and 


=0, V, = 0, 


a 


ey, 


3 
Consequently when One Ox, is applied to both sides of the equation (18) 
3 2 


and then one lets x2 = x; and x3 = x, all but one of the terms on the right 
hand side drop out. This leads to the equation 


Vil {x1}, x4, Xa) = Vo({x1}, Xn) X1y X4, 


This argument generalizes, hence © 
4. 


(19) {x1}, {x2}, ee -) 
= Vo( {x1} ™, {x2}, x1, eee, Xe, X2, 


where x; occurs p; times in the brackets. 
7. Consider now determinants similar to the Vandermondians except 
that a column is omitted. The following notation will be used. 


Form 


n—1 


(—1)*y' Vi, e(x1, Xn) + (—1)*y"** Vo, Xn) 


= Vi(y, Xn) Vo(y, M1, Xn) Ly, Xa 


Now if p:(x1, ---, Xn) is the elementary symmetric function of degree ¢ in 
Xa, ***, 


On the other hand 


Ly, X1, Xn = yxy + 


| ! | 
or- 
he | 
d, 
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nt | 
he | 
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where 


Hence 
Vi(y, Xn) = Vo(y, xn) Ly, M1, Sa 


n+k min (y, n) 


y=0 B=max (0,7—k) 
Thus for0 <t<n—1, 


(21) +, %n) = Vole, Xn)(—1)! n- 


8=max (0,t—k) 


Since the previous discussion of multiple frequencies also applies here 
we have 
Lemma 5. 


(22) Vee({ar} ™, 
t 
= Vo({xi}, ---, (— 1)? 
8=max (0,t—k) 
8. The error in using the formula h[aoF(t) + --+ + daiF(t — (m — 1)h)] 
tth 
for f F(t) dt can be explicitly computed for the case F(#) of the form e™. 
t 
Suppose that this error e(¢, e™) is in the form he(Ah)e™. Then, 
(23) + aye) 4+ + — he(dh)e* 


A corresponding equation holds for \ = »;, 7 = 1, ---,, with e(v;h) = 0. 
Cancel the e*t from equation (23) and e’* from the corresponding equations 
in the v;. Then 

Qo + + + — = — 1) 

Qo + + --- + = (hn)(e* — 1) 


It is convenient to introduce x = e~™, x; = e~** and f(x) = (1 — x)(xInx)—. 
The above system of equations is then solved for ¢(Ah). 


1, f(x1) 


Let x -—1=y, x; —1 = y;. Since Vo(x, ---, xn) depends only on 
— Xj, Volx1, = «++, Yn). The x; in the first » columns of 


ait---+an=m 
t 
4 


tions 
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the second determinant may be replaced by y;, by simple manipulations 
on determinants. 


Let g(y) = y((1 + y) In (1 + y))“. Thus f(x) = — g(y) and 
1, * g(yn) 
g(y) is analytic in y at y = 0; hence 
g(y) = Ao + Ary + Avy’ + 


When this expression for g(y) and the analogous expression for g(y,;) are 
inserted above, the elementary properties of determinants permit us to 
cancel out the first m terms in the g series. If the distributive property of 
the last column is used, one obtains for |y| < 1 and |y;| < 1, 


= - yi) (An + AnsrSi nti + AnsoSe nui + 
when we use LEMMA 3 and let 


The A, are readily calculated from the expression for g, 


y 


If the integral is expanded by the binomial theorem and integrated term 
by term, 
( 1)” 1 


A,= r(r + 1)---(r + [m — 1]) dr. 


n! 


The values up to Ay have been given in section 3. 
THEOREM I. Let do, ai, «++, @n—1 be chosen so that 
t+h 
when F(t) = e’* fori = 1, ---,m. Let e(t, F) denote 


Then for F = e and e(t, F) = he(Ah)e™, 


e(Ah) = (—1)"* (y — ys) An + nga + + 
where y = e-™ — 1, y; = e~"* — 1, provided |y| < 1 and |y;| < 1, and 


: 
—t,n- 
n> 
1)h) 
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The case of multiple frequencies is obtained by direct substitution. 
(Closed integration theorem.) Let bo, bi, ---, be chosen 
so that 


hl bo F(t + h)+ bi F(t) + b.F(t— h)+ (n— 2)h)]= F(r) dr 
when F(t) = e’t fori = 1,2, ---,n. Let e.(t, F) denote 


+ F(t — (n — 2)h)] F(r) dr. 
Then for F = e and e.(t, F) = he-(Ahje™, 


ec(Ah) = TT (y — + + 
where the y and y; are as defined above, and 
By = Ax + Ax-t. 


The discussion for closed integration is analogous to that given above 
for open integration, except that one has an additional factor e in the 
coefficients of the unknowns in the equivalent of equations (24). It is 
desirable to divide by this factor and solve for e~e,(Ahz). As a result of this 
division f. becomes (1 — x)/In x and g,(y) is y/In (1 + y). The result €.(Az) 
then is obtained by replacing the open g(y) = go(y) by g-(y) and multi- 
plying by e”. Thus the A; are to be replaced by the B; which are obtained 
from 

Bo + By + + = y/In (1 + 9). 


Inasmuch as (1 + y)go(y) = ge(y) it follows that 
B, = An+ An = (—1)" (r+ 


n! 


Thus B,, is in sign to A,, since 0 > > 


|B 


~n-1 
9. An expression is needed to represent the error e(t, 0) = sh which 
occurs when the function o(¢) is openly integrated. The following process 
will not be carried through for a closed integration procedure. Let oo = a(t), 
o; = o(t — h), etc. Then 


t+h 


This equation and the set (4’) permit one to solve for a; and s. 


Oo, ***y On-1y S 


1, Xn, 1, fly) 


e(t, F) = h[boF(t + h) + F(t) + F(t — h) + 
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To introduce the y;, let 70 = oo, m = 61 — o0, 72 = 62 — 20; + oo and 
in general 


k 
(26) m = t+ (5 + + (—1)*o». 


Note that if the o alternate in sign and are all equal in absolute value then 
m = 2*¢. With this notation, 


10, M1, ***s Nn—1s —s” 


i, ** g(yn) 
If the determinant is expanded in terms of the first row, s” appears 
| multiplied by Vo(y1, ---,¥n). However, 72 is multiplied by a determinant 
without an @ power column of the y; but with the g column. Expanding 
the g; of this column in a y series the a column cannot be cancelled as before, 
but instead one obtains the Vandermondian itself. Thus 


M0, Mis ***s 0 


n— 


s=s" — Am + 
k=0 


1, Yar D Aga! 
j=n 


We may write the error then as 


n—1 
(27) s=s"— Arm (—1)*+ o’. 
k=0 
Let 
n—1 
(28) D= ). 
k=B B 
Since 
k 
as 
B=0 B 
n—1 
(27’) s=s" — ¥ Dgog + (—1)"*' 0’. 
B=0 
Now 
M1, °° Nn—ty 0 


j=n 


n—l 
Vin, 
t=0 


A; = (— 1)'n: (- j—t—n+m, -) 
t=0 


j=n m=max (0,t+n—j) 


x 
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when Lemna 5 is used. Then 


n—-1 t 
t=1 m=0 


j=n+t—m 
Define | 
(29) G, = > 
k=0 
Then 
n—1 t 
t=1 m=0 
Thus if 
t 
(30) Tat = 
m=0 
we have 
n—1 
mT nt 
t=0 


since = It was shown that 
t 
B=0 B 
hence 
nm—1 ¢ 


B 


THEOREM II. Jf 
= pm oe 


Ge = D 
k=0 
Dy ts the elementary symmetric function of degree pin yi, Yn; 


Tat 
m=0 
n—1 t 


th 


w 


= 
3 
tl 
T 
Pp 
F n—1 n—1 T 
» 
b=0 
: 
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then the error in integrating an arbitrary function o is given by 


1 t+h n—1 
t k=0 


n—1 
t=0 b=0 


where 


and 


o’ can now be overestimated as follows: 
Let Y, 0 < Y <1 be a bound on (|y|, |y;|). |Az| < 1. Hence, from 
the definition of G,, and the result of Lemma 2 


This is a bound independent val r. As r increases, the bound may be im- 
proved by the factor | A,|. 
Now |T,,:| is also bounded. 


: ") y= 
< 
This may be inserted into the equation for H,,s, 


since (5) 


Thus 
n 
Now 


n—1 
= 
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Let ¢ = max |o;|, then 
(¥ + — 1) 


Coro.iary. The error e(t, «) of integrating a function with a maximum 
value of o is bounded by 


lo’| < 


Y+i1)"-1 
(31) < oh| 2 + — 1)| 
Project Cyclone P. Brock 
Reeves Instrument Corp., N. Y. F. J. Murray 


Columbia University 


[Editorial Note: A second part of this paper, to appear in the next issue of MTAC, will 
contain illustrative numerical examples of the application of the above ideas.] 


1W. E. MILne, “The remainder in linear methods of approximation,” NBS, Jn. Research, 
v. 43, 1949, p. 501-511. This gives a more general approach to step errors of integration 
formulas based upon approximation by sets of functions. 

2 R. E. GREENWooD, ‘‘Numerical integration of linear sums of exponential functions,” 
Ann. Math. Stat., v. 20, 1949, p. 608-611. 

3P. Brock & F. J. Murray, “Planning and error analysis for the numerical solution 
of a test system of differential equations on the IBM sequence calculator,” Cyclone Report, 
Reeves Instrument Corp., New York 28. See also F. J. Murray, ‘Planning and error 
considerations for the numerical solution of a system of differential equations on a sequence 
calculator,” MTAC, v. 4, p. 133-144. 

4F. J. Murray, “Linear equation solvers,” Quart. Appl. Math., v. 7, 1948, p. 263-274. 

5L. H. Tuomas of the Watson Scientific Computing boratory indicated ‘this formula 
for A, to the authors. He also indicated that the A, are equal in absolute value to the 
coefficients of the Adams-Bashforth method of step by step numerical integration. 

6 W. FELLER, Probability Theory. New York, 1950, v. 1, p. 52. 

7G. Birkuorr & S. MacLane, A Survey of Modern Algebra. ‘New York, 1948, p. 424. 


A Note on the Inversion of Matrices by 
Random Walks 


In a recent note, ForsYTHE & LEIBLER' described a method (first sug- 
gested by J. v. NEUMANN and S. M. ULam) for the inversion of certain 
types of matrices by a “Monte Carlo’ sampling procedure. The authors 
explain their scheme in terms of drawing balls from an urn, but the pro- 
cedure might, of course, be just as well described as a random walk. 


A boundary value problem involving a difference equation in a bounded | 


domain is equivalent to a system of linear algebraic equations in as many 
unknowns as there are lattice points in the domain. It is therefore to be 
expected that the sampling methods for the solution of such difference 
equations as explained in Curtiss? and Wasow’ are closely related to the 
method of Forsythe & Le‘bler.! 

In order to study this relation we rephrase the latter method in the 
language of random walks. We consider a set of m points P;, ---, Pm and 
introduce a moving particle which, starting from P;, jumps from point to 
point in such a way that the probability of going from P, to P, in one jump 


is P. Also at each point P, there is a probability p, = 1 — > DP» of the 


u=l 


random walk ending there. 


a 
a 
: é 


um 
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Furthermore, the moving particle possesses a variable ‘“‘mass’’ V which, 
at a step from, say, P, to P, is multiplied by a factor »,,. The initial mass 
at P; is one. Our procedure consists in estimating the expected value of the 
random variable G;; defined as follows. The moving point is known to start 
from P;: 

if the walk ends at k j 
ta Vp;", if the walk ends at 7 


Observe that G;; is defined only for points where p; ¥ 0. 
Let A be the matrix with elements 


Gig = 
(no summation is implied), denote by b;; the elements of the matrix 
B=I-A 


and by £§,; the elements of B-'. Then the following theorem is proved in 
Forsythe & Leibler.! 

Theorem: If and only if the eigenvalues of the matrix {\a;;|} are less than 
one in absolute value, then the mathematical expectation E[G,;| exists and 


E [Gi] = Bis. 


Thus, an experimental estimate of the expectation of G,; yields a numerical 
value for one element of the inverse matrix of B. 

The procedures followed in Curtiss? and Wasow’ to find Green’s function 
for difference equations, when properly worded, are special cases of a scheme 
for matrix inversion which differs from the random walk just described only 
in that the random variable G;; is replaced by the random variable M,; 
which is by definition equal to the total amount of mass carried through the 
point P; on the several visits in the course of a random walk starting from P;. 
(If the point stays put at P;, this is to be counted as a new visit.) It is 
very easy to show directly that 


= E[M,;], when p; # 0. 


For let V be the mass of the particle when it passes through P; at the end 
of a path from P; to P; whose probabil-ty of being taken is g, then 


= = Lav = E[My], 


the summation being extended over all possible paths connecting P; and P;. 
Observe that the word ‘‘path”’ is used here in a somewhat generalized sense, 
referring not to a geometric configuration but to an ordered sequence of 
points P, beginning with P; and ending with P;. 

If p; = 0, then G;; is not defined. But E[M;;] exists and is equal to B;;, 
as can be proved exactly as in Forsythe & Leibler.' This is one advantage 
of using the random variable M,; instead of G;;. Apart from this remark it 
is not easy to decide in advance which of the two methods is preferable in 
a given problem, since this requires a comparison of the variances. The 
decision is made more complicated by the fact that the factorization of the 
given number a;; into the product ,,0,; is, to a large extent, arbitrary. If 
the original problem is a boundary value problem for a difference equation, 
M;; is, to say the least, the intuitively more natural random variable, since 
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one would like to associate the end of a random walk with the first crossing 
of the boundary of the given domain. With this interpretation, p; is zero for 
all points from which the boundary cannot be reached in one step, and the 
random variable G,; is unsuitable. 

In the light of the present discussion some of the proofs in Wasow® can be 
modified—but not substantially shortened—by making use of the criterion 
for existence of E [G;;], and hence of E[M;;], which is stated in the theorem 
of this section. 

In the especially simple case that all v;; are equal to one, some additional 
information concerning the respective advantages of using the random 
variables M;; or Gi; is contained in the following theorem, valid in this 
special case. Let v; be the probability that a particle known to start from P; will 
never return to P;. Then 
(1) o[M,;] < o [Gj] 
af and only if 
Pp a < 2 v; ° 

In order to prove this inequality, we denote by \,; the probability of 
going from P; to P; without passing through P; on the way; i.e., Ax; is 
the total probability associated with all paths connecting P; and P,, all 
intermediate points being different from P;. In our special case the random 
variable M;,; is the number N of visits at P; during a random walk starting 
at P;. Then, for k > 1, Pr {N = k} = djAj*v;. A short calculation yields 


1+ 


— 


(2) E[(M,?] = E[N*] = = 


But since the total probability of the walk ending eventually is 1 we have 


(3) Aj 
so that (2) simplifies into 
1+ Ay 


Next, let X be the random variable which is 1 if the last point of the 
random walk is P;, and zero otherwise. Then 
Pr {X = 1} = + + + = — 
Hence 
E[X*] = Pr {X? = 1} = Pr {X = 1} = Agpj/(1 — ry). 


Since, in our special case, 
Gi = X/ Pi, 

it follows that 

(5) E(G?] = 


— Aw) 


ive 


the 
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Now, as the means of M;; and G;; are the same, the inequality (1) is equiva- 
lent to E[M;,?7] < E[G;?], i.e., by formulas (4) and (5) to 
1 — ry 
Application of formula (3) transforms this into p; < »;/(2 — »;), which 
proves our statement. 

The practical value of this theorem is limited, because »; is, in general, 
not known. But it shows that, at least in this special case, the answer to 
the question which method is preferable for the calculation of the element 
B,; does not depend on the subscript 7. It also confirms the intuitively 
plausible conjecture that M;; is the better random variable to use whenever 
p; is comparatively small. 


W. R. Wasow 
NBSINA 


1G. E. ForsytHe & R. A. LerBier, “Matrix inversion by a Monte Carlo method,” 
MTAC, v. 4, 1950, p. 127-129. 

2 J. H. Curtiss, “Sampling methods applied to differential and difference equations,” 
Seminar on Scientific Computation, Proc., Nov. 1949, p. 87-109. IBM, New York. 

?W. Wasow, “Random walks and the eigenvalues of elliptic difference equations,” 
NBS, Jn. Research, v. 46, 1951, p. 65-73. 


RECENT MATHEMATICAL TABLES 


961([B, E].—R. C. SpPENcER & G. E. REyNo.ps, A Table of Normalized Para- 
bolic Coordinates and Arc Lengths. Air Force Cambridge Research Lab- 
oratories, publication E 4083, Cambridge 1951, 9 p., 21.9 X 27.3 cm. 


The functions 
y = 4x2, s = Harcsiah x + x(1 + x*)!) 


are given to 5D for x = 0(.01)2. These give the ordinates and arc distances 
of points on the “normalized” parabola y = $x*. From these tables two 
graphs are derived showing s as a function of x and s, x as functions of y. 
From the latter graph the excess of arc length over abscissa may be read 
off to facilitate the laying out of a parabolic antenna. 

D. H. L: 


962[F].—J. P. Kuiix, L. Potettr & R. J. Porter, Liste des Nombres 
Premiers du Onziéme Million (plus précisément de 10006741 @ 10999997). 
[Amsterdam 1951, published by N. G. W. H. Beeger, Nicolaas Witsen- 
kade 10, Amsterdam]. ii + 25 p., 20.3 X 27.6 cm. Price 3 Dutch florins. 


This table is the result of collating three of the four existing manuscript 
factor tables of the eleventh million.! The arrangement is essentially that 
of LEHMER’s? list of primes of which this table is a natural extension. Thus 
the rank of a prime occupying page P, column C, and line L is given by 


2500P + 100C + L + 662400. 


The number of primes in this million is 61938. 
Much of the credit for the successful completion of this table goes to 
BEEGER and GLODEN who were responsible for collating and reconciling the 
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manuscripts and for the preparation of the final clear handwritten manu- 
script from which the table was offset. 
D. 8. L. 


- 1See MTAC, v. 4, p. 121, for a complete description of this project. 
2D. N. Lenmer, List of Prime Numbers from 1 to 10006721. Washington 1914. 


963[H].—B. M. SuumfAcsxil, Tablisty dlia resheniia kubicheskikh uravnenii 
metodom osnov [Tables for the solution of cubic equations by the method 
of bases]. Moscow and Leningrad, 1950, 135 p., 12.9 & 20.1 cm. 3.90 
rubles. 


In a previous paper! the author discussed properties of the solutions of 


(1) 2¥+Az—A=0, 
and their application to the solution of the more general trinomial equation 
(2) m+ py" +q=0. 


In particular when m = 3 and nm = 1, then N = 3 and (2) is transformed 
into (1) by letting 
(3) y=-—(q/p)z and A = p*/¢. 


The resulting cubic equation (1) has three real roots 2: < 2. < z3 when 
A < —6.75; as A decreases to — ©, these vary monotonically over the 
ranges 2; = (—3, —©), 2. = (1.5,1), 23 = (1.5, ©). When A > —6.75, 
(1) has one real root, which is the continuation of z,, and two complex roots 


(4) 2,3 = — (2:/2) + ina; 


2, increases from —3 to 1 and a from 0 to © as A becomes infinite. 

The booklet under review has a dozen pages of explanation and examples 
followed by four tables which contain for the most part values of A (and a 
in the complex case) for equally spaced 2, (and z. when real). There are 
about ten pages with values of 2; (or 22.) and a (when the roots are complex) 
for equally spaced A. Non-trivial A are given for all functions. 

Table I contains A and a for 2; = .0020(.0001).010(.001).969 and 2; and 
afor A = 30(1)300. Table II gives A and a for 2; = —.0020(—.0001) —.010- 
(—.001)—3. Table IIIa has A for 2} = —3(—.001)—7 and z, for A = —40- 
(—1)—30. Table IIIb shows A for 2, = 1.5(—.001)1.030 and z for 


A = —37(—1)—300. For |A| > 300 the following formulas are suggested 


without comment by the author (p. 9): 
A>300, a=(4+2)/(A+3), = —(2:/2) + V—A — 
A < —300, = (A+ 2)/(A+3), 21 = —(@/2) + V—A — 


The bulk of the values of A, z:, and z, are computed to 5S while most 
of those of a are given to 4S. On p. 6 the author writes, ‘“The tables give 
accuracy to four decimal places for real roots and three decimal places for 
complex roots. The 4th and 5th places for the first case and the 4th place 
for the second case can be obtained by linear interpolation. In must cases 
interpolation will give 5 correct places of a real root. Pages where linear 
interpolation gives only four places (but not five) of a real root have an 
asterisk in the upper left corner.” 
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83 


NocGrapy’ has tabulated to 6D the functions A = 2/(1 — 2) for 
a = —3(.001)1 and A = 2z,*/(1 — 2) for 2, = 1(.001)1.5. [Actually he uses 
the variable Z = —z, or —z:.] Consequently, his values of A are computed 
to more decimal places than those on p. 19-26 (.400 < 2 < .971), p. 45-70 
(—1 > 2, > — 3), and p. 126-131 (1.500 > z > 1.030) of the present volume. 

If equation (2) (with m = 3, m = 1) is to be transformed into a one- 
parameter equation, the substitutions (3) leading to the form (1) (with 
N = 3) require only rational calculations with p and g in contradistinction 
to the irrational computations needed to put (2) into the more popular forms 


|j+2z| =a and +2 = 3pz. 


ZAVROTSKY® has tabulated explicitly both real and complex solutions of 
(2) to 5D for p = —100(1)100, g = 0(1)100. 
R. R. REYNOLDs 
NBSINA 


1B. M. SaumiAcsktl, “‘Ischislenie shoirosh,” Mat. Sbornik, v. 40, p. 394-409, 1933. 

2H. A. NoGrapy, A "New Method for the Solution of Cubic Equations. Ann Arbor, 1936 
{MTAC, v. 1, p. 441-442]. 

3A. ZavRrotsKy, Tablas para la Resolucién de las Ecuaciones Ctibicas. Caracas, 1945 
{MTAC, v. 2, p. 28-29]. 


964[K].—NiLs BLomovist, ‘“‘Some tests based on dichotomization,” Annals 
Math. Stat., v. 22, 1951, p. 362-371. 


Consider a random sample of m vectors for an m-dimensional population 
with unknown distribution. The only allowable values for a component of 
a vector are 0 and 1. Thus the vectors can be represented in the form 
(Xi, Xam) (¢ = 1, where each x,; is either 0 or 1. Let 


a = pj = n/nj, (j = 1, ---, m) 


x4, = 1,-+-,m). 
j=1 
The y; are considered to be random variables but the ; (consequently, 
also the p;) are considered to be fixed. Let 


bi S= — PY, 
i= i= 

while all the x;; are considered to be statistically independent. Also let 
(u,v) = (n = u,m =v). Table I presents values of P(S = So) for: (4, 3), 
So = 1(2)9; (4, 4), So = 0(2)10, 16; (4, 5), So = 1(2)13, 17, 25; (4, 6), 
So = 0(2)20, 26, 36; (4, 7), So = 1(2)13, 17, (2)21, 25(2)29, 37; (4, 8), 
So = 0(2)26, 30(2)40, 50; (6, 3), Sp = 1.5(2)9.5, 13.5; (6, 4), Sp = 0(2)18, 24; 
(6, 5), So = 1.5(2)25.5, 29.5; (8, 3), So = 0(2)14, 18; (8, 4), So = 0(2)26; 
(10, 4), Sp = 0(2)30; (10, 3), Sp = 2.5(2)18.5; (12, 3), So = 3(2)23; (14, 3), 
So = 3.5(2)25.5; (16, 3), Sp = 0(2)28. Values given have either 3D or 4D. 
Number of significant figures =3 (except for 1.000); 4D used when number 
of significant figures =2. 


J. E. Watsn 


U. S. Naval Ordnance Test Station 
China Lake, California 
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965[K, P].—_E. BRockMEYER, Sandsynlighedsregningens Anvendelse i Tele- 
fonien [Applications of the Calculus of Probability to Telephony]. 
Copenhagen, 1949, Copenhagen Telephone Company. 54 p., 16.8 X 25.0 
cm. 


This work contains a table of solutions y to 2D of the equation 
(1) BY y*/k! = y*/n! 
k=0 


for B = .001(.001).005, .01, .02, .05 and for m = 1(1)260. The table is 
intended to indicate the traffic intensity y which n circuits will handle with 
the probability B of loss of a call. The relation (1) is known as Erlang’s 
formula [MTAC, v. 3, p. 98-99]. This table appears also in RMT 966. 

D. H. L. 


966[K, P].—E. Brockmeyer, H. L. HAtstrgm & ARNE JENSEN, “The life 
and works of A. K. Erlang,”’ Akad. Tekn. Videns., Copenhagen, Trans., 
1948, no. 2, 277 p. 


There are small tables related to the Poisson distribution included in 
this work. Among those of some general use are 


(a) 6D values of m*e~"/x! for x = 0(1)14 

and m = —2(.1)0. This table is on p. 137. 
(b) 6D values of (—m)*e"/x! for x = 0(1)20, m = 2.1(.1)40n p. 195-196. 
(c) 6D values of the real and imaginary parts of the solutions of 

z exp (z) = a exp (—a) exp (27ik/40) 

for a = 0(.1)1, & = 0(1)20 (p. 197-198). 

(d) BROCKMEYER’S table, described in RMT 965. 
D. H. L. 


967[K].—M. B. Cook, “‘Bi-variate k-statistics and cumulants of their joint 
sampling distribution,” Biometrika, v. 38, 1951, p. 179-195. : 


This paper has for its purpose the derivation of ‘‘(a) formulae giving 
bivariate population moment-coefficients in terms of cumulants, and cumu- 
lants in terms of moment-coefficients up to the sixth order; (b) formulae 
for the cumulants of the joint distribution of the simpler bivariate k-statis- 
tics.” In pursuance of (a) the author indicates three methods of obtaining 
the results but he actually uses Kendall’s operational method! to put 
on record the bivariate moments about zero and about the means through 
order 6 in terms of cumulants and conversely the bivariate cumulants 
through order 6 in terms of moments about the means. In the second part 
of the paper he turns to the calculation of the cumulants of the joint distri- 
bution of bivariate k-statistics.? This is carried out by an extension of the 
operator method used in the first part. He lists results through weight 9. The 
results, being largely new, were checked by calculation in two different ways. 


J. J. Livers 
The Boeing Aircraft Company 
Seattle, Wash. 


1M. G. KENDALL, “The derivation of multivariate sampling formulae from univariate 
formulae by symbolic operation,” Annals Eugenics, 10, 1940, p. 392-402. 
2R. A. FisHer, ‘Moments and product moments of sampling distributions,” London 
Math. Soc., Proc., s. 2, v. 30, 1929, p. 199-238. 
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968[K].—S. T. Davin, M. G. Kenpat & A. Stuart, “Some questions of 
distribution in the theory of rank correlation,” Biometrika, v. 38, 1951, 
p. 131-140. 


Let (Xi, Xo, ---, Xn) be a random ordering of the integers (1, 2, ---, 2), 
and let S(d*) = é (X; — i). The authors give (Tables 1, 2) the exact 


probability distribution of S(d*) for n = 9, 10. [The distribution for 
nm = 1(1)8 has been published earlier, e.g. by KENDALL.'] These tables permit 
exact tests of the significance of the Spearman rank correlation coefficient 
for samples of size n. By obtaining algebraic expressions for the moments 
and cumulants through order 8, the authors provide approximate distribu- 
tions for S(d?) which are sufficiently accurate for practical purposes when 
n = 10. 

J. L. Hopces, Jr. 
University of Chicago 
Chicago, Ill. 

1M. G. KENDALL, The Advanced Theory of Statistics. v. 1. London, 1948, p. 396. 


969[K].—(a) J. Dursin & G. S. Watson, “Testing for serial correlation 
in least squares regression. II,’’ Biometrika, v. 38, 1951, p. 159-178. 
(b) G. S. Watson & J. Durstn, “Exact tests of serial correlation using 
noncircular statistics,’ Annals Math. Stat., v. 22, 1951, p. 446-451. 


Tables in (a) facilitate the test for serial correlation of error terms, e, 
of a regression model y = 61x; + Box2 + --- + Bex, + € where y and the x; 
are observed and the x; regarded as “‘fixed variables.” Let z denote the 
residual from regression and Az denote the successive first differences of z. 
The statistic d = }-(Az)?/>°2* is used for the test. Define k’ = k — 1 and 
set x, = 1. Upper and lower bounds (dy and dz) are tabulated for d, where 
P(d > daz) = a for a = .05, .025, .01 for sample size m = 15(1)40(5)100, 
for k’ = 1(1)5 and assuming normality of residuals. The statistic d is 
always between 0 and 4 and tabled values are given to 2D. The use of 
these tables causes the test to be inconclusive if d, < d < dy and approxi- 
mate procedures are recommended for this case. Examples of the use of the 
above tables are given for one independent and two dependent observations, 
two-way classification and polynomial regression. 

Two small tables appear in (b). The first is 95th percentile values for 
the distribution of c, where 


for n = 2m and where xms1 = 0 for n = 2m + 1 (effectively making m = 2m) 


given to 3D for m = 10(2)22. The second table gives 95th percentile values 
for the distribution of d; where 


for n = 2m and samples of 2m + 1 are reduced to samples of size 2m by 
omitting the middle observation. Values to 3S are given for m = 12(2)30. 
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Much of this table appears in Table 7 in (a) which also considers upper and 
lower bounds for the 95th percentile of a comparable statistic for least 
squares regression for k’ = 1(1)5 and m = 13(2)23. 

W. J. Dixon 
University of Oregon 
Eugene, Ore. 


970[K].—A. K. GayeEn, ‘‘The frequency distribution of the product-moment 
correlation coefficient in random samples of any size drawn from non- 
normal universes,’ Biometrika, v. 38, 1951, p. 219-247. 


The author’s summary includes the following description of part of the 
content of the paper and of Table 1 of the paper. ‘“The mathematical form 
of the distribution of r in non-normal samples is obtained, the parent popu- 
lation being specified by the bivariate Edgeworth surface including terms as 
far as those in Ago”, AgoAe1, +, AosArz aNd ---. The tables of the ordi- 
nates of the corrective functions due to fourth-order and the square and 
product of third-order semi-invariants . . . are obtained for VN = 11, p = 0.0, 
and 0.8, and N = 21, p = 0.8.” 

These corrective functions provide multipliers for certain functions of 
the semi-invariants \,;. In terms of formulas, Table 1 gives values to 5D of 


dt 
(cosh t — 


V(r, p) = (1 — p*)#4—-D (1 — 72) f 


N-1 
p) = + 1) apt v(r,e) for «= 1,2, and 
N-2 ai 


p) 12N(N + 1)(N + av’ p) for «= 1, 2,3. 


¥(r, p) is the ordinate of the distribution of the correlation coefficient 
in samples from a bivariate normal population and its values are taken 
from F. N. David’s table.1 The derivatives of y with respect to p are ex- 
pressible simply in terms of y's. Values are given for N = 11, p = 0.6, 
r = —1.00 (.05) 1.00; VN = 11, p = 0.8, r = —.75 (.05) .60 (.025) 1.000; and 
N = 21, p = 0.8, r = —.25 (.05) .60 (.025) 1.000. In the last two cases the 
entries are zero to 5D for smaller values of r than those given. 

The author finds the semi-invariants of several bivariant samples which 
exist in the literature and discusses, in part with the aid of tables and graphs, 
the effect on the distribution of r of taking into account semi-invariants of 
order greater than two. The effects of non-normality in a bivariate popula- 
tion of the distribution of z = 3 log. [(1 + r)(1 — r)~] is also discussed. 


K. J. ARNOLD 
University of Wisconsin 


Madison, Wis. 


1F, N. Davin, Tables of the ordinates and probability integral of the distribution of the 
correlation®coefficient in small samples. London, 1938. 


971[K].—D. G. C. Gronow, “Test for the significance of the difference 
between means in two normal populations having unequal variances,” 
Biometrika, v. 38, 1951, p. 252-256. 
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The significance of the difference is to be tested by one of the two’ 


criteria 
1\f (m — 1)s:? + (m — 


v= — + =)" 


Ne 


where s? = }- (xu — &)?/(m: — 1), t = 1, 2. Only two-tailed tests are con- 
i=1 


sidered. One sample of m, observations is taken from a normal population 
with mean £ = x, and variance x, and one sample of m, observations is 
taken from a normal population with mean & = «;’ and variance x,’. The 
tables in the paper give probabilities that each of u and » falls beyond the 
5% level and that each falls beyond the 1% level of ¢ (two-sided) for 
m, + m, — 2 degrees of freedom, when m = m, = 10 and when m = 15, 
m2 = 5 for various combinations of 6/ Vig and ko’ /ke. 6 = & — &. Proba- 
bilities are given for m, = m2 = 10 (in which case u = 0), 5/Vk2 = 0(.5)3, 
ka! /x = 1(.5)3.0, and for m = 15, = 5, = 0(.5)3.0, = 1/3, 
1/2, 1, 2, 3. 4D are given for 6 = 0, 3D otherwise. 

As the author states, “this problem has already been solved in principle 
by Hswu,} but his solution is not in a form which readily allows of numerical 
calculation.”” The author has therefore approximated the first four moments 
of the distribution of each of his statistics for each case considered, fitted 
“the appropriate Pearson or Gram-Charlier Type A curve,” and used this 
approximation to the distribution of the statistics to obtain the tabled 
probabilities. Checking by using Hsu’s results in a few cases in which they 
can be simplified, no serious discrepancy was found although it appears 
likely to the reviewer that the 7 in 0.276 for m = m, = 10, 1% level, 
8/Vk2 = 1, x2’/xe = 1 is a misprint. When 6/x. appears in a heading or in 
the text it should be 5/Vx:. 

K. J. ARNOLD 
University of Wisconsin 
Madison, Wis. 


1P. L. Hsu, “Contribution to the theory of ‘Student's’ é-test as applied to the problem 
of two samples,” Stat. Res. Mem., v. 2, 1938, p. 1-24. 


972[K].—E. S. KEeEptne, “A significance test for exponential regression,” 
Annals Math. Stat., v. 22, 1951, p. 180-198. 


The author considers the two regression equations 
(1) Y = be and (2) Y=a+ be 


under the assumptions that the values of x are in arithmetic progression 
and that the standard deviation of the observed y is constant for all x. The 
null hypothesis being tested is that be?* = 0 for all x while the alternative 
hypotheses are specified by b ¥ 0, p # — ~. In Table III the critical values 
of R, the sample correlation coefficient (uncorrected for means) between the 
observed values of x and y in the case of (1), are given for samples of 3 (1) 10, 
12, 15 at the 5%, 1%, .1% and .01% levels of significance to 3D at least. 
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In Table IV for the range » = 4 (1) 12, and the same significance levels as 
in Table III, the author provides the critical values of R to 3D at least, for 
(2), where R in this case is computed in the usual way. 


L. A. AROIAN 
Hughes Aircraft Co. 
Culver City, California 


973[K].—L. R. Packer, ‘‘The distribution of the sum of m rectangular 
variates I,’’ Inst. Actuaries Students’ Soc., Ju., v. 10, 1951, p. 52-61. 


The author rederives, using difference methods, the well-known distri- 
bution function for the sum of 2 randomly drawn values from the continuous 
rectangular universe on the interval (0,1). The cumulative distribution 
function, F(x,m) is tabulated to 6D for m = 1(1)12 and x = .5(.5)6. As 
an application he gives a table of the exact probabilities that the error in 
the sum of 12 numbers, each rounded off to the nearest unit, will lie in the 
interval (—x, x) for x = .5(.5)6. 


974[K].—E. S. Pearson & MAXINE MERRINGTON, “Tables of the 5% and 
0.5% points of Pearson curves (with argument §; and 8.) expressed in 
standard measure,’’ Biometrika, v. 38, 1951, p. 4-10. 


As an alternative to using the normal distribution as an approximation 
to the exact distribution of a statistic, the authors propose that one might 
use the distribution of the Pearson system corresponding to the moment 
ratios, 

Bi = Bo = pape? 


of the statistic. This proposal could be adopted if one can compute the first 
four sampling moments of the distribution, and if tables of areas for the 
corresponding Pearson curve are available. 

As a start towards the satisfying of the latter requirement, the authors 
have tabulated the values of 


u = (Xa — where f f(x) dx =a 


for a = 0.005, 0.05, 0.95, and 0.995, in which yp, o, 6; and ® are the moment 
ratios of the statistic whose distribution one is interested in, and f(x) is the 
solution of the corresponding Pearson differential equation, 


dy/dx = — (x + (Go + + 


The values of (x — u)o— are tabulated to 2D for 6; = 0, .01(.02).05(.05).2 
(.1)1, and for B. = 1.8(.2)5. 

By using these tables one can approximate the 5% and 0.5% critical 
points of a statistic’s distribution. Tables which will yield additional critical 
points are in process of being computed. The appropriateness of using the 
Pearson system for such approximations is briefly discussed but no criterion 
which would help one determine when he should use this method is given. 


C. F. Kossack 
Purdue University 


Lafayette, Ind. 
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975[K].—K. C. S. Prixat, “On the distribution of an analogue of Student’s 
t,”” Annals Math. Stat., v. 22, 1951, p. 469-472. 


Let & be the mean and w the sample range in a sample of m from a nor- 
mal population with mean a and variance unity ; further define G = fond. 
The author lists in Table I for m = 3 (1) 8 (2) 14 (3) 20, (1) the confidence 
interval for a based on the G test, (2) the confidence interval for a based on 
Student’s ¢ test, and (3) the efficiency of the G test, defined as the ratio, 
(1)/(2), all to 4D. For this range the efficiency of the G test is always greater 
than 96%. The levels of significance, a, are .1, .05, .01 and .001 for nm = 3; 
.1, .05, and .01 for m = 4 and 5; and .05 and .01 for the other values of n. 


L. A. AROIAN 


Hughes Aircraft Co. 
Culver City, California 


976[K].—M. S. Rarr, “The distribution of blocks in an uncongested stream 
of automobile traffic,” Am. Stat. Assn., Jn., v. 46, 1951, p. 114-123. 


Cars passing through a fixed location on a highway are observed; the 
assumption is made that the number k of car passages during a time interval 
of length ¢ has the Poisson distribution e~¥*(N#)*(k!)—. A car arriving at 
this location on a side road is assumed not to enter the highway until there 
is a “gap” of duration =Z between consecutive passages of cars on the 
highway. Intervals of time containing no gap of duration =L are called 
“‘blocks”’ of traffic, while each gap of dugation =Z is called an ‘“‘antiblock.” 
Various probability distributions relating instants of arrival of cars on the 
side road (assumed to be uniformly distributed), blocks, antiblocks and 
their durations are derived. The main result of the paper is the probability 
distribution of the ‘‘waiting time’”’ of a car arriving on the side road, from 
the moment of arrival to the earliest antiblock. The cumulative distribution 
function of the waiting time is 


r=0 


where h is the greatest integer contained in ¢/L. This distribution depends 
on{two parameters NL and t/L. A table of F(t) is presented for NL = 0(.25) 
2(.5)6, t/L = 0(.5)3(1)7(2)11 and selected larger values until F(t) reaches 
.99. Properties of F(t) are discussed. An empirical set of data is compared 
with a theoretical distribution and, on inspection, shows good agreement. 
Z. W. BIRNBAUM 
University of Washington & Stanford University 
Seattle, Washington & Palo Alto, California 


977[K].—P. R. Riper, “The distribution of the range in samples from a 
discrete rectangular population,” Amer. Stat. Assn., Jn., v. 46, 1951, 
p. 375-378. 


The author considers the infinite statistical universe in which the integers 
0, 1,2, +--+, NW — 1, occur with equal frequency and derives the probability 
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that in a random sample of m from it the range will be R. The values of 
this probability are tabulated to4D for N = 10,” = 2(1)10, and R = 0(1)9. 
For comparison the like table is given for samples of m from the continuous 
rectangular distribution on the interval (0, 10). 

<<. 


978[K].—S. S. SHRIKHANDE, “Designs for two-way elimination of hetero- 
geneity,’’ Annals Math. Stat., v. 22, 1951, p. 235-247. 


Designs are presented using ten or less replications for the following 
numbers of treatments and treatments per block. 


Treatments per block Number of treatments 


6, 10, 13, 15*, 25 
8, 10, 15 


CONA US W 
= 
an 


10 31 


Methods are also presented to construct designs for s* treatments in blocks 
of size s. 


R. L. ANDERSON 
North Carolina State College 


Raleigh, North Carolina 
* Using either 5, 6 or 7 replications. 


979[K].—C. H. Sprincer, ‘‘A method for determining the most economic 
position of a process mean,” Industrial Quality Control, v. 8, 1951, p. 
36-39. 


The discussion and tables given here deal with the problem of proper 
location of the process mean for the case in which the loss due to producing 
a piece above the upper specification is not equal to the loss of producing 
a piece below the lower specification. 

Let X be a random variable (quality measure) with (process) mean X’ 
and (process) standard deviation o’. Let ¢ be the transformed random 
variable with mean zero and standard deviation one. Let U be the upper 
and L the lower specification, Cy the cost of rejecting a piece X > U and 
Cx the cost of rejecting a piece X < L. We have for C(X’) the total ex- 
pected unit cost of rejected product, 


C(X’) = Ci iy f@dt+ of, f(b) dt. 


The first condition for C(X’) a minimum ie 
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which, together with the second minimizing condition, gives conditions on 
the density function f(t) which are satisfied by (among others) the Gaussian 
and Pearson Type III considered by the author. For the former (1) re- 
duces to 


= L+ Wo’, where W C1/Cu, Cu < Cu, 


o’ 


L 


W and W’ being equal for reciprocal cost ratios. Springer tables W for 
fifteen selected values of Cz/Cy (eight values, including unity, and their 


= U— W'e', where W’ = -( ) in C1/Co, Cu > Co, 


reciprocals) in the range .005 to 200, and for = _— = 2.5(.5)6.0. 


Similar tables are given for the Pearson Type III function. Letting 


K 
n (o’)? 
where m is the sample size, Springer gives tables of W for K = —0.3, 


—0.7, —1.1, 0.3, 0.7, 1.1, for eleven selected values of Cz/Cy (six values, 
including unity, and their reciprocals) in the range .01 to 100, and for 
— 2.5(.5)5.5. 


Haroip A. FREEMAN 


Massachusetts Institute of Technology 
Cambridge, Mass. 


980[K].—D. R. Wuitney, “A bivariate extension of the U statistic,’’ Annals 
Math. Stat., v. 22, 1951, p. 274-282. 


Let x1, +++, Xt} ***, Zn be independent variables having 
(by hypothesis) the same continuous cumulative distribution. Let U be 
the number of times a y precedes an x and V be the number of times a z 
precedes an x when the observations are arranged in order of size. Table I 
gives the number of sequences of the x, y and 2’s having specified U and 
V(U = 0(1)18, V = 0(1)9) for the special case ] = 6, m = n = 3. Table 2 
gives 1000 times the cumulative distribution of U and V (and also a bi- 
variate normal approximation for this) for the same special case and range 
of values for U and V. 

The author states that tables of the joint distribution of U and V have 
been prepared for all cases where 1 + m +n ¢ 15. 


FRANK J. MASSEY 


University of Oregon 
Eugene, Oregon 
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981[L].—MiLton “Table of the integral, f dn,” 
0 
Math. Physics, v. 30, 1951, p, 162-163. 


The integral is tabulated for x = 0(.1)2.5 to 8D, with second central 
differences. For x < 1.5 the computation was carried out from the power 
series, for x > 1.5, by numerical integration. In both cases 10D were ob- 
tained and independent checks were made. Linear interpolation is said to 
be generally good to 4D. 

A. E. 


982[L].—I. BLocu, M. H. Hutt, Jr., A. A. BRoyLes, W. G. Bouricrus, 
B. E. FREEMAN, & G. Breit, ‘‘Coulomb functions for reactions of 
protons and alpha-particles with lighter nuclei,’ Rev. of Modern Phys., 
v. 23, 1951, p. 147-182. 


In many calculations in theoretical physics, the differential equation 


occurs. Its two most important solutions are the ‘‘regular Coulomb func- 
tion,’’ F,, and the “irregular Coulomb function,’’ Gz, which are determined 
by their asymptotic behavior for large (positive) p, 


Fy, ~ sin (9 — — 2p + oz), 
Gr ~ cos (p — 4La — 2p + oz) 
where 
o, = arg T(L+1+%). 


There are also several auxiliary functions. Constants Cz, Dz are defined by 
Co = — 1) L2L+1)Cr = + (224+ = 1; 
and functions Az, gz, Ox, by 


Fr, = At sin = Crp™6,, = 
Gr = Arcos = Drp-* Oz, Gz’ = 0; *. 


The aim of the numerical tables presented here is to provide values, 
with an accuracy of one percent, required in the problems mentioned in 
the title. 

There are 39 tables, pages 155-182. Tables I to V: Fy to Fy. Tables VI 
to X: &) to &. Tables XI to XV: auxiliary functions for the computation 

- of Sp to & in ranges in which Tables VI to X cannot be interpolated to 1 
percent accuracy. Tables XVI to XIX: o* to #;*. Tables XX to XXIV: 
to Tables XXV to XXIX: to Tables XXX to 
XXXIV: go to gy. Tables XXXV to XXXIX: Ap to Ay. 

login = —.8(.1).6 in the tables, and 0 < p < 6.2 at varying intervals. 
Not all values of logio and of p occur in all tables. The values are given 
to 4D, and each table is accompanied by a note on accuracy and inter- 
polation. 

References to previously available tables are given on page 152. Methods 
for computation of the functions tabulated here are described in another 
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paper’ by the same authors, and on pages 153-154 of the present paper three 
numerical examples illustrate the use of the tables. 


A. E. 
1]. Brocu, M. H. Hutt, Jr., A. A. Broytes, W. G. Bouricius, B. E. Freeman, & 


G. Breit, ‘Methods of calculation of radial wave functions and new tables of Coulomb 
functions,” Phys. Rev., v. 80, 1950, p. 553-560. 


983(L].—_R. J. BuEHLER & J. O. HirscHFELDER, “Bipolar expansion of 
coulombic potentials,” Physical Rev., v. 83, 1951, p. 628-633. 


P, and P, being two points of three-dimensional space at distance ry 
from one another, let P; be referred to a system (r1, 61, $1) of spherical polar 
coordinates with pole O,, and P2 to a system (re, 62, ¢2) with pole Oz, it being 
understood that the polar axis is common to the two systems (and passes 
through O, and O,): also R is the distance 0,0>. 

In the expansion 


1 
72; R)Pm, (cos 6:)Px,(cos 62) cos m(¢2 — ¢1) 
ni1,n2,m 
the coefficients B can be represented in simple closed forms if either 
< Ror |r; — r2| > R. In the remaining case, — < 1, 
the B assume a more complicated form, 


2(ny+nq+1) 
1% 


In table I, the coefficients D and A are given exactly for m, m,. = 0(1)3 and 
for all appropriate values of m,i, 7. The results should be useful in the 
computation of the interaction of two overlapping Coulomb fields of 
potentials. 


A. E. 


084[L].—_F. R. Erskine Cross.ey, “A hyperelliptic function as a non- 
linear oscillation,” Jn. Math. Phys., v. 30, 1951, p. 214-225. 


Tables I-V give, for g = 10°, 45°, 90°, 135°, 170°, ¢ = 10°(10°)90°, 
and e? = .10, .25, .5, .75, .90, 1, values to 4D of ¢ in radians, sin ¢, and of 


1 — e(1 — 2k? sin? 
f [ 1 — sin? Ja. k= 


Appended are values, in degrees, minutes, and seconds, of ¢ defined by 


sin = sin sin 
2 2 
There are also some graphs connected with the integral. 


A. E. 


985[L].—W. R. Dean, “Slow motion of viscous liquid in a semi-infinite 
channel,” Cambridge Phil. Soc. Proc., v. 47, 1951, p. 127-141. 
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Among the several short tables in this paper one finds 4D tables for 
u = .9(—.1)—.9 of z 


1 2 1 —i1i— 


2u(1 + u?) w(1 + u)4 u(i + u?) + u)4 


(itu)? 16u(1 + 


1-—u 


= ip — 
where / = log (1 — u), = (2%) 


A. E. 


986[L].—Otto EMERSLEBEN, “‘Numerische Werte des Fehlerintegrals fiir 
Vnx,” Zeit. angew. Math. Mech., v. 31, 1951, p. 393-394. 


There is a table on p. 393 of 15D values of 


f exp (—#) dt 


for x = (xn), m = 1(1)11. 

R. C. ARCHIBALD 
Brown University 
Providence, R. I. 


987[L].—A. N. Gorpon, “The field induced by an oscillating magnetic 
dipole outside a semi-infinite conductor,’ Quart. Jn. Mech. Appl. Math., 
v. 4, 1951, 106-115. 
Table I (p. 110) gives exact values of r~* and 4D values of the real and 
imaginary parts of 


t3) 


for r = .2(.2)1(.5)2. Here b = it. 
A. E. 


988[L].—T. S. Kuun, “A convenient general solution of the confluent hyper- 
geometric equation, analytic and numerical development,” Quart. Appl. 
Math., v. 9, 1951, p. 1-16. 


The confluent hypergeometric equation is written in the form 


[- 2 + 1) 
dr + r r 


and z = (8r)!. The author puts 


| Un = 0, 


Un) U,© (z) 
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d 
and gives explicit expressions of U,(z) and 4z ie U,(z) for 1 = 0,1 and 


k = 0(1)7 in terms of Bessel functions of order 0 and 1: the coefficients of 

the Bessel functions are polynomials in z. Two particular solutions arise if 

the Bessel function is J or Y, and for these particular solutions the explicit 

expressions were evaluated, and are tabulated to 4D, for z = 3.5(.5)7.5. 
A. E. 


989(L].—R. O. Gumprecat & C. M. S.iepcevicu, Tables of light-scattering 
functions for spherical particles. xvi + 574 p., $6.50. 
Tables of Riccati Bessel functions for large arguments and orders. xvi + 260 
p-, $3.50. 
Tables of functions of first and second partial derivatives of Legendre 
polynomials. xii + 310 p., $3.50. University of Michigan Engineering 
Research Institute, Special Publications. Ann Arbor, 1951. 15 K 23 cm. 


The diffraction of electromagnetic waves on a sphere has been investi- 
gated by several authors, most completely by G. Mreg,' and the theory is 
presented in standard text-books.? The infinite series occurring in Mie’s 
theory converge the slower the larger the radius of the sphere (in comparison 
to the wave-length of the incident light-wave). The NBSCL tables* provide 
data for spheres whose radius does not exceed the wave-length of the incident 
light. From the practical point of view much larger droplets are also of 
interest, and the present tables supply data for spheres whose circumference 
is up to 400 times the wave-length. On account of the slow convergence of 
the series, it is estimated by the authors that the computation of the tables 
with a standard desk-computer would have taken 25 man-years. Using 
IBM equipment this period could have been reduced to 3 man-years. 
Actually, the computation was carried out on the ENIAC, and less than 
two weeks’ time was required once the details of programming had been 
worked out. 

The notation adopted is that of the NBSCL tables, and is explained in 
greater detail in MTAC, v. 3, 1949, p. 484. The relevant parameters are 
a, m, y and the following notations are used: 


aP,,(x) @P,,(x) 
dx’ 


where P,(x) is the Legendre polynomial, J the Bessel function of the first 
kind, and S, and C, are called Riccati-Bessel functions. There are also two 
complex quantities A, and P, which are compounded of Riccati functions 
and their derivatives, the variables being a and 8. [This P, should not be 
confused with the Legendre polynomial in the definition of x, and ,’. ] The 
quantities of physical interest are the intensity functions 


B=ma, x=cosy, = 


= {Anta + — (1 — x*)x,']} 


n=1 


2 
h= 


n=1L 
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and the scattering function 
A,|? + |P,|? 
K(m, 0) = 3 2 Gat 
[In the definition of the intensity functions the square brackets in place 
of | | are presumably a misprint. ] 

In the first volume Part I gives values (mostly to 4S) of the intensity 
functions for y = 90°, and of the scattering function, for m = 1.2, 1.33, 1.4, 
1.44, 1.5, 1.6 and a = 1(1)6(2)10(5)100(10)200(50)400. Values of the inten- 
sity functions for other angles y are not tabulated. They may be computed, 
however, from the values of An, Pr, tn, etc. Part II of the first volume 
gives 6D values of the real and imaginary parts of A, and P, for the same 
a and m as above, and for a suitable range of . The practical difficulties 
in computing the series for larger a (and 8) are illustrated by the following 
examples indicating various ranges of m appropriate to the different values 
of a, and m = 1.2: 


a n 
1 1(1)4 
10 1(1)17 
100 1(1)114 
400 1(1)421. 


In the second volume under review those Riccati-Bessel functions are 
tabulated which were required for the computation of the A, and P,. 
Values (mostly to 6S) are given of S,(z), S,’(z), Ca(z), Cn’(z), in Table I for 
z = a with a and a ranging as above, and in Table II for z = ma, with 
m, a, and m ranging as above, except that those values which have already 
appeared in Table I are omitted. Certain other omissions, especially for 
m = 1.33 and m = 1.6, also occur. Judging from the manner in which the 
infinite series converge, the authors believe that the tabulated values are 
accurate to within one unit of the last decimal place given. 

The last of the three volumes under review contains material for the 
computation of intensity functions. Specifically, it gives tables of 7, and of 
xm, — (1 — x*)z,’ to 9 or less S for m = 1(1)420 and y = 10°(10°)170°(1°) 
180°. The tabulated functions are stated to be accurate to within one unit 
of the fifth significant figure, except when fewer than 5S are given in which 
case one unit of the last place is the maximum error. 

The first of the three volumes has a foreword by P. DEeByE. Each 
volume has an introduction giving definitions of the tabulated quantities, 
a brief description of the computation, notes on accuracy and interpolation, 
and a request that users communicate to the authors any errors or omissions 
noticed. In the reviewer’s copy there is one correction in ink: i, for a = 35, 
m = 1.33, y = 90°, should read 9.926 instead of 14.07. 

The tables should be very useful in practical computations of scattering 
from droplets. On the whole they are well arranged and well produced, and 
perhaps it is ungracious to complain about details. There are a few blem- 
ishes, though, most of which could have been avoided. In the second volume 
the table headings were stencilled, and it is difficult to distinguish between 
C and c. Also, the tables in this volume have been reduced considerably, 
perhaps excessively, in the photographic process. In the third volume the 
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table headings were typed, and the same symbol was used for the variable x 
and the multiplication cross X. 
A. E. 
1Gustav Mie, “Beitrage zur Optik triiber Medien,” Amnnalen der Physik, s. 4, v. 25, 
1908, p. 377-445. 
2]. A. STRATTON, Electromagnetic Theory. New York, 1941. 


’ NBSCL, Tables of Scattering Functions for Spherical Particles. Washington, D. C., 
1948. See MTAC, v. 3, 1949, p. 483-484. 


990[L].—F. W. J. Over, “A further method for the evaluation of zeros 
of Bessel functions and some new asymptotic expansions for zeros of 
functions of large order,” Cambridge Phil. Soc., Proc., v. 47, 1951, 
p. 699-712. 


In a previous paper,’ the author described a method for determining 
zeros of Bessel functions J,(x) and Y,(x). The method is independent of 
tabular values and zeros are obtained by numerical integration of a third 
order non-linear differential equation. The order of the functions is fixed 
and integration is performed with respect to a variable ¢. The numerical 
procedure was not satisfactory for the early zeros and an alternative ap- 
proach was developed. Here the technique amounts to the computation 
and inverse interpolation of the phase of the Hankel function. 

The first part of the present paper gives a new method for rapidly 
determining the early zeros. The idea is to treat ¢ as fixed and regard the 
order of the Bessel function m as a continuous variable. If p(n, ¢) is a zero 
of the real cylinder function C,(x), then the starting point is an equation 
of WaTson,? p. 508, 


(1) = F(p,n) = 2p f ” Ko(2p sinh u)e-®™ du. 


The asymptotic expansion of (1) is determined for large p and m. With 
= n/p, the result is 


Fp, 2) = (— 


(2) = (— f sech*"*! uam(A sech u) du. 


Values of a, are given by Watson,? p. 316, for m = 0, 1, 2, 3. In the appli- 
cations 0 < A < 1, and it is advantageous to tabulate f,,(A) over this range 
so that interpolation is easily accomplished. For this purpose, the author 
defines 


(3) In(A) = f (cosh u + du. 


Values of J; and J: are easily obtained and J,, is determined recursively. 
The f,, are then expressed in terms of (3). The author derives a formula so 
that values of the derivative of a function at its own zero can be found by 
quadrature. The formula involves an asymptotic expansion whose coeffi- 
cients depend on gn(A) = (—1)""f,,’(A) which are written in terms of (3). 
In a similar fashion expansions are developed so that zeros of C,’(x) and 
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stationary values of C,(x) can be obtained by numerical integration. Here 
again the necessary coefficients depend on (3). 

To facilitate the calculations, 10D values of J,,(A) have been computed 
for m = 1(1)13, A = 0(.05)1. These values were then used to compute 
Sm(A), Zm(A) for m = 0(1)4 and A = 0(.01)1. These tables are not given in 
the present paper. The procedures of the paper are being used to evaluate 
the zeros of J, Y and their derivatives as part of the program of the Royal 
Society Mathematical Tables Committee (see MTAC, v. 3, p. 340). The 
reviewer anticipates that tables of f,,(A), etc., will be published simulta- 
neously with zeros of J, Y, etc. To illustrate techniques of paper, extracts 
of tables (1 and 2) are given showing the numerical integration process for 
computing the first zero of J,(x) to 10D over the range » = 6(.5)20. 

In the second part of the paper, asymptotic expansions are given for 
the direct evaluation of the following: 


(a) positive zeros p of the real cylinder function C,(x) 
(b) values of C,’(p) 

(c) positive zeros ¢ of C,’(x) 

(d) values of C,(c) 


The first few coefficients needed in (a) to (d) are tabulated in Tables 3 to 5, 
respectively. Values are given corresponding to the first five zeros of J,(x), 
Y,(x), etc. Asymptotic expansions are given for item (b) and also its recip- 
rocal. It is for the latter that the coefficients are given in Table 4. Coefficients 
for the former can be easily evaluated in terms of Table 3 and a portion of 
Table 4. Accuracy obtained using tabulated coefficients is judged by ex- 
amples. If s = 1, m = 20, seven-figure accuracy can be obtained from items 
(a) and (c); six-figure accuracy, in associated values of items (b) and (d). 
If s = 5, = 20, the corresponding results are five and four-figure accuracy. 
Accuracy increases with m and decreases with s. 

The basic coefficients needed for the above tabulations depend on tables 
of the Airy integral.* The present tabulations are virtually all new and 
supersede the earlier work of MEIssEL‘ and Arrey® whose results for the 
most part are in error. Airey’s results are quoted in JAHNKE & EmpE.*® 

YUDELL L. LUKE 
Midwest Research Institute 
Kansas City, Missouri 


1F, W. J. Otver, “A new method for the evaluation of zeros of Bessel Functions and 

of 46, 1950, reo of second-order differential equations,” Cambridge Phil. Soc., Proc., v. 
GN. Watson, A Treatise on the Theory of Bessel Functions. Cambridge, 1945. 

s BAASMIC, Pt. v. B. The Airy Integral. Cambridge, 1946. 

‘E. MEISSEL, “Beitrag zur Theorie der Bessel’schen Functionen,” Astr. Nach., v. 128, 
1891, cols. 435-438. 

5J. R. Arrey, “The numerical calculation of the roots of the Bessel Function J,(x) 

and its first derivative J»'(x),”” Phil. Mag., s. 6, v. 34, 1917, p. 189-195. 

SE. JAHNKE & F. Empe, Tables of Functions. New York, Dover, 1945, p. 143. 


991[L].—K. M. SrEGEL, D. M. Brown, H. E. Hunter, H. A. ALPERIN, & 
C. W. QUILLEN, The Zeros of the Associated Legendre Functions P,.™(u’) 
of Non-Integral Degree, University of Michigan, Engineering Research 
Institute, Report No. UMM-82, Ann Arbor, Mich., 1951, vii + 20 p. 
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In a previous discussion of a paper by CaRRUs and TREUVENFELS (CT), 
a difference test indicated that some of the early zeros of the associated 
Legendre function P,'(cos 6) = 0 as a function of m were incorrect (MTAC, 
v. 5, p. 152-153). The investigation of the present article also reveals 
some errors. The authors give an alternative proof of an equation due to 
MACDONALD? for determining the early zeros of P,™(cos 0) = 0 where @ is 
near 7. For m = 1, 6 = 165°, this formula gives 1.035 as an approximation 
to the first zero. Employing power series, it is shown that the first zero 
must be between 1.0316 and 1.0321. The value reported by CT is 1.053 
and so is in error. Application of the Macdonald formula shows that for 
165° < 6 < 180°, the corresponding values of m decrease with increasing 0. 
For m = 1, @ = 170°, the first zero given by CT is 1.05 and thus is also 
incorrect. Numerical analysis of early zeros for values of @ other than those 
cited above is not given, but sufficient evidence now exists to show that 
the CT tables should be used with caution. 


L. LUKE 
Midwest Research Institute 
Kansas City, Missouri 
1P. A. Carrus & C. G. TREUVENFELS, ‘‘Tables of roots of incomplete integrals of asso- 
ciated Legendre functions of fractional orders,”’ Jn. Math. Phys., v. 29, 1951, p. 282-299 
MTA v. 5, p. 152- 153}. 


. MACDONALD, ‘ ‘Zeros of the spherical harmonic P,”(u) considered as a function 
of n,’ ’ London Math. Soc., Proc., s. 1, v. 31, 1900, p. 264-278. 


MATHEMATICAL TABLES—ERRATA 


In this issue references to Errata have been made in RMT’s 989, 990, 
and 991. 


203.—AKADEMIiA Nauk, SSSR. Tabliisy znachenit Funktsit Besselia ot 
mnimogo Argumenta. |MTAC, v. 5, p. 151-152.] 


p. x Function For Read 
10 444 Ail 2367 2267 
19 .899 Ai 667 657 
42 2.031 A, 593738 493738 
42 2.032 A, 481922 381922 
106 5.237 1H 153939 132939 
114 5.650 AA, 788 782 
115 5.700 % 5.605 5.700 
118 5.815 AM, 476 469 
118 5.816 AH, 449 456 
163 8.061 A, 949 849 
166 8.235 AM, 001 81991 
195 9.654 1H 276029 276022 
205 074 AK, 82290 81290 
206 139 AKo 76 876 
220 815 AKo 693 683 
220 848 AKo 645 685 


230 1.312 Ko 380745 381745 
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p. x Function For Read 
238 1.719 J; 20163 30163 
249 2.295 J; 889777 869777 
251 2.369 J4 942091 941991 
253 2.491 J; 318266 308266 
254 2.538 J; 480753 490753 
254 2.539 J; 506447 516447 
257 2.663 J; 864568 884568 
260 2.803 Ky 926830 926820 
264 3.044 AK, 565 505 
264 3.047 AK, 394 334 
269 3.262 J4 739537 739587 
269 3.293 J4 271376 271326 
276 3.629 AKy 122 128 
293 4.480 AJ; 779 784 
294 4.502 AKy 588 582 
304 5.029 AKy 145 140 
374 8.533 Ko 300300 300250 
375 8.576 K, 118097 118197 
392 9.423 AKo 283 263 


Besides the above errata, 112 errors of less than five units in the last place 
were noticed. For further errata see MTAC, v. 5, p. 152. 


S. A. JOFFE 
515 W. 110 Street 


New York 25 


204.—BAASMTC, Mathematical Tables, v. 1. Cambridge, 1946, 2nd ed. 


Table II—Circular Functions 
page 7, x = 48.6 
for .0945447099 79701 read .0945447098 79701 


STANLEY H. CoHN 
Indiana University 
Bloomington 


205.—N. W. McLacuian & P. HumBeErt, Formulaire pour le Calcul Sym- 
bolique. Mémorial des Sciences Mathématiques, fasc. 100, 1941. 


.4, formula 2; upper limit of second integral: for t read wt*/2. 

. 5, formula 11; omit the index ¢ on the left hand side. 

. 7, formula 3; same correction to second integral as on p. 4. 

. 14, formula 2; omit b+ on the right hand side. 

. 14, last formula; for 1 on the r.h.s. read p.. 

32, formula 10; to the I.h.s. add (1/2) log [(¢ — 6)/(t + 6) ]Jo(ay). 
34, formulae 3, 4; to the l.h.s. add +(i/7) log [(¢ — 6)/(t + 6) ]Jo(ay), 
respectively. 

34, formulae 5, 6; to the l.h.s. add + (it/m) log [(¢ — b)/(t + 6) JJo(ay), 
respectively. 

. 34, delete formulae 7, 8. 

. 55, formula 1; for m above >> read m. After the formula add m = 4n, 
n even; m = 3(n — 1), n odd. 
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p. 58, formula 1; should be on p. 19 in § 3, since |sin ¢| is not discon- 
tinuous. 
p. 61, line below Régles; for fo(x) read fe(y), and for $2(p) read ¢2(q). 
N. W. McLacuLan 
Vizard & Co. 
51 Lincoln’s Inn Fields 
London, W.C. 2 


206.—N. W. McLacuian, P. HumBert & L. Supplement au Formu- 
laire pour le Calcul Symbolique, Mémorial des Sciences Mathématiques, 
fasc. 113, 1950. 


4, formula 4; upper limit of second integral; for x read wx? /2. 
6, third formula from bottom; for Vix read ivix. 

7, last formula; same correction as on p. 4. 

8, formula 4; forO < x < reaadO <x < 

. 10, line 19; see correction to p. 14 in Fascicule 100. 

. 11, line 11; see corrections to p. 34 in Fascicule 100. 

. 19, last line but one; for K read k. 


ov 
. 29, formulae 3, 4; for —~—— = read sin (aV3t/2). 


. 30, formula 3; for 1 in i. cee read p. 

. 44, formulae 9, 10; for read 

. 46, formula 2; for e~ read e~*. 

. 47, formula 4; the argument of the second gamma function should be 
(—u —» + 3). 


D 


N. W. McLacuLan 


207.—E. OBerG & F. D. Jones, Machinery’s Handbook. 14th edition, New 
York, 1950 (and earlier editions). 
Table of Gear Ratios and Decimal Equivalents, p. 708-711 (different 
page number in earlier editions) 
Insert in appropriate positions 


Decimal Equivalent Gear Ratio 
.6944 25/36 
.7956 34/45 
Tuos. H. O’ BEIRNE 
Barr & Stroud 
Glasgow, W. 3 


208.—K. Pearson, Tables of the Incomplete Beta- Function. Cambridge, 1934. 


On pages 2 and 3, line 1, p = g = 3 in the value of the complete Beta 
function 
for 3.14159245 read 3.14159265 


H. W. Norton 


University of Illinois 
Urbana, Illinois 


102 AUTOMATIC COMPUTING MACHINERY 


209.—J. V. Uspensky, Introduction to Mathematical Probability, 1937. 
On page 407, Table of the Probability Integral 
for $(z) = .499997 read .500000 


CHARLES T. JOHNSON 
5852 Adelaide Avenue 


San Diego, California 


UNPUBLISHED MATHEMATICAL TABLES 


In this issue an Unpublished Manuscript Table is referred to in RMT 
990. 


144[F].—A. GLopEN & J. BoNnNEAU, Factorization of N*+ 1 for isolated 
values of N betweer 30000 and 40000. One page typewritten manuscript. 
Deposited in UMT FILe. 


The table contains 88 factorizations, all complete. No primes are given. 
[For previous tables of this kind see MTAC, v. 2, p. 211, 252, 300; v. 3, 
p. 21, 118-119, 486; v. 4, p. 24; v. 5, p. 133-134.] 


145[D, F].—D. H. Leumer, Table of Cyclotomic Cosines. Ten manuscript 
pages tabulated from punched cards. On deposit in the UMT FILE. 
Also available on punched cards. 


The table gives 20D values of 
2cos2rk/p for k = 1(1)(p — 1)/2 
for every odd prime p < 100. There are 517 values in all. Thus the table 
gives twice the real parts of the p-th roots of unity. 


146[F, L].—D. H. Lenmer, Table of Kloosterman Sums. Twenty manuscript 
pages tabulated from punched cards. On deposit in the UMT FILE. 
Also available on punched cards. 


The table gives 19D values of 
p—l 
Sp (k) = exp {2mi(kn + 7)/p} (nn = 1(mod p)) 


for k = 1(1)p — 1 and for every odd prime p < 100. The table was com- 
puted from UMT 145, and contains 1034 entries. These sums appear in 
Fourier coefficients of many elliptic modular functions. 


AUTOMATIC COMPUTING MACHINERY 
Edited by the Staff of the Machine Development Laboratory of the National Bureau 
of Standards. Correspondence regarding the Section should be directed to Dr. E. W. 
Cannon, 415 South Building, National Bureau of Standards, Washington 25, D. C. 
TECHNICAL DEVELOPMENTS 
THE SERIAL-MEMORY DIGITAL DIFFERENTIAL ANALYZER 


Introduction. In January, 1950, the first model of a digital differential 
analyzer became a working reality. This machine was entirely contained 
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in a volume of 5} cubic feet. The compactness, simplicity of construction, 
and digital accuracy of this computer by far overshadowed any differential 
analyzer built prior to this time. 

This first machine possessed 22 integrators, each of 22 binary digit 
capacity. 

The key to the practicality of this relatively midget computer is its 
serial memory which enables one network of computing circuitry, which 
operates on one binary digit at a time, to be shared for all computations. 
This computing center consists of 44 vacuum tubes and 700 diodes. 

Although any type of serial memory can be used for this type of machine, 
the most practical memory to date is a magnetic drum because of its com- 
pactness and reliability. 

It is the purpose of this paper to give a brief explanation of the logical 
embodiment of this type of computer and a description of some of the 
machines most recently designed. 

General Description. The Digital Differential Analyzer can be consid- 
ered in four basic sections: The control unit, the memory, the computing 
center, and the power supply. 

The control unit provides for filling the memory with information, 
starting and stopping computation and for displaying memory information 
by means of a readout indicator. This unit also controls power. The unit 
can be built as a part of the computer or it can be provided as a separate 
unit for remote operation. 

The memory is a rotating drum coated with a black magnetic oxide. 
With the exception of a permanently recorded clock channel, all memory 
channels are dynamically recorded between a read and a record head covering 
a sector of the drum. Appropriate amplifiers and shapers are associated 
with the read and record heads to develop proper signals. 

Since the computer is a serial device, information, as it is read from 
the memory, passes through the computing center a pulse (binary digit) 
at a time. Here it is altered as demanded by the calculations performed and 
then re-recorded in the memory. The computing center is a maze of gates, 
inverters and flip-flops which are closely interrelated so that all computation 
is performed in synchronism with the memory clock. 

Integrators. This computer, like all differential analyzers, is made up of 
fundamental integrating units, known as integrators. An individual inte- 
grator accepts two streams of digital information as a series of input pulses 
and releases, as an output, pulses related to the inputs by integration. Since 
the inputs and outputs of integrators are in the same physical form (pulses), 
a problem is solved by making proper electrical connections between 
integrators. 

Consider the two registers and additive-subtractive transfer arrange- 
ment shown in figure 1. This is in essence a digital integrator in which 
the dy and dx of the integrator are represented by pulse rates. Likewise 
the output of the integrator, dz, is in the same form. If Y is added to 
R every time a pulse occurs on the dx line and Y is created by accumu- 
oo dy increments, this device operates such that the pulse rate out 


= K( dy a) = K where is the rate of pulses accumulated 
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; 
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in Y and “ is the rate of Y to R addition. K is dependent on the number 


system used and the numerical capacity of R. [For more details see MTA C, 
v. 6, p. 41-49.] 

The Memory. Recirculation or delay-type channels are used in this type 
of computer for the entire memory. A permanently recorded clock channel, 
which is closed upon itself around the drum, is used for synchronism through- 
out the computer and memory. (All reading and recording of memory 
signals and computer operations are kept in step by pulse times defined by 
these clock signals.) 

Each incoming signal from the computer unit to the memory is first 
amplified and clipped in an amplifier whose plate is clamped to allow a 
limited swing. This insures that signals of uniform amplitude and wave- 
form are applied to the grids of the recorder tubes. The ‘‘non-return-to-zero”’ 
method of recording is employed, whereby a record head is used_to reverse 
a saturation pattern previously laid down by an erase head. 


R 


Q 


| 
| 
dt | 
| 
| 
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ot 


Fic. 1. Transfer device which performs digital integration. 


The recorded pattern is picked up by the playback head after,an appro- 
priate delay represented by the fraction of a revolution between the record 
and playback heads. After detection, the signal is amplified and sent to the 
computer section. Here it is put into a phase inverter which gives two 
amplified outputs, both similar to the input signal but opposite in phase. 
Since the detected signal from a magnetic memory is the derivative of what 
was recorded, a positive pulse represents a change of 0 to 1 and a negative 
pulse from 1 to 0. The negative pulses from the two voltage sources of the 
phase inverter correspond now to these two memory changes. By passing 
the phase inverter signals through clippers arranged to pass only negative 
signals, the signals which denote changes of memory content are indicated 
by well-defined negative pulses. These pulses are then put unto the appro- 
priate grids of a flip-flop (primary memory flip-flop) and the memory infor- 
mation is reconstructed in the same ‘‘non-return to zero” waveform which was 
originally recorded. Because of the “‘record-detect-reconstruct” process just 
described, there is no guarantee that the memory waveform will agree with 
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the clock-defined pulse times. To eliminate any relative jitter that may be 
present, the primary memory flip-flop gates clock synchronized signals to a 
second memory flip-flop which reconstructs the original recorded signal, 
delayed by a specific number of clock pulses. 

The registers used in the transfer system described in the example of 
digital integration are all segments of the magnetic drum memory. Elec- 
tronic circuitry needed is thus reduced to a minimum, since the same circuits 
of the computing section are employed for all registers on a time-sharing 
basis. The R & Y numbers are recorded on R and Y channels respectively. 
These channels are divided so that appropriate capacity is allotted to any 
one set of registers. 

The indexing of these registers, recorded sequentially in the memory, 
is performed by means of counters. Assume the computer under discussion 
is to represent m integrators (mY registers and mR registers), each of which 
has m pulses (binary digits) of capacity. By the count in an electronic 
counter of capacity , any specific pulse under consideration, of a particular 
register, is then located in time. The individual registers are similarly 
accounted for in a counter of capacity m, which keeps track of the inte- 
grators. It can be seen that the memory capacity is m X n pulses, which is 
the factor determining the spacing between read and record heads. 

Interconnections between integrators are made by using a third channel 
on the drum; this is called the Z or precessing line. The Z channel is used 
to store the dz output pulses of the integrators as they pass through the 
computing center. This is accomplished by recording the overflows of the R 
registers which occur at the Pn, or last pulse time, of an integrator. 

The Z line is much shorter than the other memory channels, its length 
being relative to the number of pulses in one integrator time. By making 
the delay of this line one pulse time longer or shorter than the number of 
pulses of one integrator, the pattern recorded will shift by one pulse for each 
successive integrator time, allowing the recording of new overflows to take 
the place of obsolete dz information. As long as there are more pulse posi- 
tions in the Z line than there are integrators, all integrator outputs can be 
stored in the Z channel and will be available to all integrators. 

Corresponding to the R and Y channels there are two address channels 
Dx and Dy. Proper coding of these channels enables the associated inte- 
grator to receive dx and dy information from other integrators. This is done 
by coincidence circuits. Since the length of the Z line and the timing and 
ordering of the information supplied to it are known, the time that specific 
dz information is available for pickup by a given integrator can be predicted. 
By placing a pulse in the proper position of an address channel, Dx or Dy, 
desired overflow information or dz outputs of any integrator contained in 
the Z channel, are available for the digital transfer process as defined by 
dx and dy inputs of one of the integrators. 

The nature of the machine prevents the direct coding of more than one 
dx pickup for any one integrator; however, more than one dy pickup is 
possible as long as the number of dy’s used for any one integrator does not 
exceed the capacity designed into the dy section of the computer. 

Block Diagram of Computer Operations. The block diagram of figure 2 
illustrates an abbreviated serial-memory (magnetic drum) digital differ- 
ential analyzer, during the computing operation. 
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Consider first the Z channel: the information in this channel is re- 
recorded through gate G1, except at Pn time when this gate closes to block 
off old information and pass a new overflow from the R+ Y adder through 
gates G2 and G3. This G2—G3 system gates an overflow, provided one occurs, 
with the exception of Pn of word W1, at which time G3 is blocked and an arti- 
ficial overflow corresponding to W1 is automatically recorded through G2 into 
the memory. In this manner a constant pulse rate dz/dt, which is the rate 
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Fic. 3. Precession line for machine of 5 integrators of 5 pulse layers each. 


of the independent variable, is generated. The Z information also enters a 
dy decoder where the Y pickup code extracts the proper pulses from the 
Z line and accumulates the proper dy increments. This information is then 
available, as a 2 dy signal, to the Y portion of an integrator. Timing is such 
that the proper ‘‘connections” reach the right integrators. 

The question now arises as to the manner in which the Z information is 
recorded to distinguish sign. Instead of a plus-zero-minus system, which 
would be impossible to record in a binary fashion in only one channel, a 
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plus-minus system is used which has the zero rate biased as alternate ones 
and zeros (i.e. if a given integrator has successive outputs of 1—0-1-0, etc., 
its output dz rate is equivalent to zero). This means that any pulse rate 
from the Z line (representing the continuous outputs of a particular inte- 
grator) that has a time average of more ones (recorded overflow) than zeros 
(absence of a recorded overflow) is positive, while the opposite (more zeros 
than ones) is a negative rate. The mechanics of this system as it is used in 
the dy decoder will be explained in a discussion of the number system. 

Corresponding to the dy decoder is a dx decoder where the Dx code 
enables an examination of the Z line to take place, resulting in a signal 
corresponding to the dz information held in the pulse position of the line 
being inspected. A dx (+) signal as an output from the decoder corresponds 
to a one in Z while dx (—) corresponds to a zero. The dx (+) signal controls 
the addition of Y to R while dx (—) controls subtraction of Y from R. 

The address lines never change content during computation and recir- 
culate as shown. These channels are used only to time proper Z line inspec- 
tion in the decoding sections for integrator interconnections. 

The first pulse in the Y section of an integrator, called the start pulse, 
is a marker which determines the size of the Y number which will be con- 
tained in the following portion of the Y number. This pulse is detected by 
the S flip-flop, which remains in the on state for the rest of the integrator 
time, turning off at Pn. The numerical portion of Y is added to = dy in 
the Y + = dy adder resulting in a new Y value, denoted as Y*. This addition 
is controlled by the S signal as indicated. Ky is a delay flip-flop holding the 
carry from one binary digit to the next. 

Gate G4 passes Y back to the memory as long as the S flip-flop is off. 
Since the first pulse occurring in Y is the start pulse, this is the only infor- 
mation other than zero passed by G4 and it is thus held unchanged in the 
Y memory during computation. Gate G5 passes the new Y, (Y*), into the 
Y channel of the memory when S is on. The Y channel then always holds 
the start pulse and the current Y + 2 dy sum of the integrators being 
used in computation. The integrators not used are left blank and recirculate 
zeros through G4. 

The R + Y* adder operates in the same manner as the Y + = dy adder. 
It has inputs of R, Y*, and a carry delay flip-flop Kr. However it does not 
have to be controlled by S since the start pulse is eliminated from Y* in 
the Y + 2 dy adder. Since dx controls the addition of Y* to R, to the 
extent that Y may either be added or subtracted, special precaution must 
be taken. This is accomplished by gates G6 and G7. If the output of the 
dx decoder signifies addition, the dx (+) signal opens G6 and allows Y* to 
be added to R; on the other hand if the dx (—) signal is present, G7 is open 
and the output of a complementing circuit is added to R (this is the same 
as subtracting Y* from R). 

The output of the R-+ Y* adder is denoted by R*. As previously 
explained, the overflow pulse of R + Y* = R* is the output of the adder 
at Pn; this is gated into the Z precession loop through G3. 

As was the case for the start pulse of the Y channel, the information 
recorded in R at Pn must remain unchanged from one memory cycle to 
the next. This effect is accomplished by G9 and G10. The pulse position 
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Pn in R is occupied by what is known as a “‘sign-reversal pulse.’”’ A one in 
this pulse position reverses the sign of the dz output. This operation can be 
explained as follows: if this pulse position stores a zero a certain pattern 
of ones and zeros as successive overflows will occur from the R + Y¥* adder 
for a specific integrator. 

Consider now the same conditions as above except with a permanent 
one existing in Pn of R. The overflow pattern results in an extra one added 
to the sum at each Pn time, thus resulting in the ones and zeros of the sum 
at this particular pulse time being interchanged, which is a reversal of sign 
of the output rate. 


NUMBER SYSTEM 
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Y Number System. The contents of Y are shown in figure 4. The first 
pulse of Y, Pi, is the first possible start pulse. The remainder of Y, to Pa, 
represents the corresponding numerical content of Y. In the case where 
the start pulse is at P;, P, to P,_; is numerical information. P, is always 
blank, for this corresponds to the overflow pulse of R. P,-: is the sign pulse; 
a one in this pulse position denotes a positive number, while a zero denotes 
a negative number. The number system, with the example of the start pulse 
contained in P,, is shown in figure 4b. This arrangement allows a Y number 
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capacity of +(2*-* — 1). In general the capacity of Y is +(2"-*-* — 1) 
where the start-pulse is placed in P;. Different arrangements of start-pulse 
positions are illustrated in figure 4. 

The dy decoder can be considered as a counter which adds the dz over- 
flow pulses in the Z line coincident with the dy code channel pulses and 
subtracts pulses whenever a dy code pulse is coincident with no dz pulse. 
Thus at the end of such an accumulation the counter holds the algebraic 
sum of the 1’s and —1’s. When this sum (called = dy) is sent to the Y + = dy 
adder, the transfer is done in a serial manner. This transfer is accomplished 
by a continuous stepping of the counter contents from the first stage Vj, 
the signal of the “‘on plate” of this flip-flop being passed into the adder 
circuit. The 2 dy number system and read-out are shown in figure 5 for a 
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computer of 3-stage accuracy. It is the size of this counter-register which 
determines the dy pickup capacity of a given design. It is obvious at this 
point that the accumulation and stepping cannot be performed simulta- 
neously. For this reason two units are used for this overall operation. The 
= dy accumulation takes place in one unit, the 2 dy counter, an integrator 
time before it is used. It is then transferred into a stepping unit, the 2 dy 
register, at the beginning of the next integrator time, and finally stepped 
into the adder circuit. Figure 6 illustrates the addition of 2 dy + Y. 
The Precession Line. To illustrate the operation of the precession line, 
figure 3 shows a line one pulse time longer than an integrator for a hypo- 
thetical machine of 5 integrators, each of 5-pulse capacity. The information 
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at the read and record heads is always the same, except when a new overflow 
is recorded at the time the new information is introduced to the Z line at 
the read head. The line is originally clear at the start of computation, 
picking up its first bit of information at the overflow time of J, (PsJ;). 
This overflow 1 progresses in the line until at P;J_ it is at position 5; at this 
time 2 is being recorded at the record head. The process builds up until 
at Psis all five overflows are recorded. At PI, (second cycle) a new 1 is 
recorded replacing the old 1, “shaded position” in actual use. At P;J, the 
old 1 is replaced by a new 2 which makes the old 2 obsolete. This process 
continues, making all previous overflow information available to integrators 
immediately following. The pattern of the Z line can now be defined for 
any pulse of any integrator. At P;J;, the information occurring at the read 
head is the last output of integrator 4 — 7 + k. The precession lines used 
in larger computers are essentially elaborations of this simple example. 

Filling and Controls. The discussion of the serial-memory digital differ- 
ential analyzer has been limited to the computing operation. Memory filling 
operations are closely related to those of computation in that some of the 
circuitry of the computer is used for filling. 

The essential controls for filling are: channel erasing, channel selection, 
integrator selection and control of information filled. The erasing control 
is built into the memory circuitry and consists of a controlled recording, 
by means of clear switches, of zeros in a channel. The channel selection 
is accomplished by a switch that controls the logical circuitry at the record 
end of the computing section. An integrator selector is wired into the 
integrator counter and matrix in such a manner that a signal is generated 
corresponding to the integrator time set by this selector. Information is 
filled into the chosen integrator and channel by keys corresponding to 
the binary digits, 1 and 0. These keys fill a corresponding one or zero into 
successive pulse positions P; to P,. 

The method of filling successive pulse positions is quite simple and 
makes use of the Z line. When the machine is idling (neither computing 
nor filling) all memory channels recirculate, including the precession (Z) 
line. Before filling, a marker pulse must be placed in the Z line by the 
operator. This pulse serves as a reference for the pulse position to be filled. 
When one of the filling buttons is pressed, a circuit, operated by the coinci- 
dence of the fill button, integrator selector signal and marker pulse, fills 
a pulse time in the selected channel with a 1 or 0 corresponding to the button 
pushed. The marker pulse then shifts by one pulse and is ready for the 
next digit to be filled. The stepping is done by allowing the Z line to precess 
once during filling operations at the word time coincident with filling. 

CRC 105 Decimal Digital Differential Analyzer. The DDDA (more 
commonly known as D#A) is the first decimal machine of this type to be 
designed. The coded-decimal system makes the simple design features of 
its predecessors applicable and in many instances introduces simplification. 
This computer will have 60 integrators, each of six-decimal place capacity, 
plus a position for indication of sign. 

External inputs to the DDDA will be stored in the precession line, and 
unnecessary coding will be eliminated thereby. The inputs will be picked 
up in the same manner as integrator outputs. Typewriter output on this 
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machine will allow the contents of up to 12 integrators to be typed out at 
an interval computed by the machine. This typewriter will eliminate much 
of the need for auxiliary equipment such as graph plotters, etc. 

The computer will be able internally to limit the maximum and minimum 
of functions, which is of great advantage in dealing with non-linearities 
occurring in electronics and servo problems. 

Means of storage of initial conditions during computation will eliminate 
the tedious refilling ordinarily necessary in repetitive problems, where 
families of solutions or experimental solutions are made. 

The serial-memory digital differential analyzer was co-invented by 
D. E. Ecxpant, R. E. Spracue, H. H. Sarxissian, C. L. Isporn, A. E. 
Wo FE, W. F. CoL.ison, and F. G. STEELE. 

Joun F. DonaN 
Computer Research Corporation 
Hawthorne, California 
DIscussIONS 


The Adventures of a Blunder 


Those who regularly code for fast electronic computers will have learned 
from bitter experience that a large fraction of the time spent in preparing 
calculations for the machine is taken up in removing the blunders that have 
been made in drawing up the programme. With the aid of common sense 
and checking subroutines the majority of mistakes are quickly found and 
rectified. Some errors, however, are sufficiently obscure to escape detection 
for a surprisingly long time. The blunder which is the excuse for this note 
falls in this class and to the writers’ knowledge was responsible for only 
one arithmetical error in its whole career! How many other errors of a 
similarly obscure character are still at large in the EDSAC library can only 
be conjectured but one at least has been reported to us by Dr. A. vAN 
WIJNGAARDEN. 

The blunder arose from the misuse of a subroutine for calculating the 
square roots of numbers in floating decimal form. Numbers in this form are 
expressed as a-10”, where p is an integer. Calculations are carried out 
using a “floating decimal accumulator,” that is, a pair of storage locations 
set aside for holding a number which is being operated on in floating decimal 
form. Before being transferred from the floating decimal accumulator to 
another part of the store, each number is standardised so that 1 < |a| < 10. 
In the floating decimal accumulator itself, |a| may lie outside this range 
(the representation being not unique). Thus in the multiplication of two 
numbers each in standard form, the two numerical parts are multiplied 
yielding a number in the range 1 < |a| < 100, and this is left in the floating 
decimal accumulator as it stands. Similarly, addition of two standard 
numbers may leave the numerical part anywhere in the range 0 < |a| < 20, 
the exponent being that of the larger component. 

A programme was prepared which involved finding the modulus of a 
complex number X + 7Y. A square root subroutine was used for finding 
the square root of the quantity X? + Y?. This square root subroutine causes 
a number held in the floating decimal accumulator in standard form to be 
replaced by its square root. However, the requirement that the argument 
should be in standard form was overlooked in the preparation of the sub- 
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routine. Strangely enough the programme passed lengthy tests successfully 
and the blunder was only discovered when the coding of a similar problem 
led to the investigation of this particular detail. Further investigations 
showed that the chance of the blunder producing an error was indeed 
small, about 1 in 256. The explanation is given below. 

First it was found that the square root subroutine would, in fact, func- 
tion correctly if |a| < 16, or if p was even and |a| < 160, and it happened 
that this restriction was met by nearly every case encountered in the pro- 
gramme. The reason for this is as follows. Let the complex number be 
X +4Y, where X is expressed as x-10* and Y as y-10°. First, X? was formed 
and transferred to the store, being thus standardised as, say, z-10°, which 
could be either x?- 10? or (x?/10)-10°*+4. Then Y? was formed in the floating 
decimal accumulator and X? added to it, and the square root subroutine 
was Called in. At this stage we can distinguish two cases: 


(i) s < 2r. The exponent in the floating decimal accumulator is 27, and 
X? + Y? appears as (y? + z-10*-*")-10". The exponent is even and the 
numerical part does not exceed 110; hence the subroutine works correctly; 

(ii) s > 2r. The exponent is s and the numerical part is (y*- 10”-* + 2)10*. 
If s is odd, s = 2g + 1 and the numerical part may be as high as 16 only 
if g =r, in which case the restriction may be violated. The subroutine 
works correctly unless g = r, or x? + y? > 160. The likelihood of this 
violation can be roughly estimated by assuming a logarithmic probability 
distribution of the modulus and a rectangular distribution of the argument. 
With this assumption the probability of an error is .0039. The error is thus 
rare enough to escape notice in all but exhaustive tests. 


University Mathematical Laboratory 
Cambridge, England 
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1. ANon., Digital Computer Newsletter, ONR Mathematical Sciences Divi- 
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The contents are as follows: 


. The OARAC 

. The CADAC 

. The SEAC 

The SWAC 

Whirlwind I 

Moore School Automatic Computer (MSAC) 
The Burroughs Laboratory Computer 

. Naval Proving Ground Calculators 
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Special Purpose Computer 
1. MADDIDA 
Data Processing and Conversion Equipment 


1. The Charactron 

2. Wallind-Pierce Equipment 

3. Aurex Magnetic Wire Storage 

4. CRC Ferro-Resonant Flip-Flop ~ 


2. ANON., “New instruments,”’ Rev. Sci. Instruments, v. 22, July 1951, p. 
544-545. 


High-speed analog to digital converter announced by Arthur D. Little, 
Inc., Cambridge 42, Mass., takes an input voltage varying from 0 to 50 
volts, samples it up to 50,000 times a second, and converts the readings to 
7-digit binary numbers expressed as pulses on seven output leads. 

Analog-digital converter announced by Glenisco, Inc., 2233 Federal 
Avenue, Los Angeles 64, Calif., counts shaft rotations in the decimal system 
to a least count of 0.1 revolution by means of solenoid actuated quantizers. 


R. D. ELBourRN 
NBS Electronic Computers 


3. Stuart R. BRINKLEY, JR., ‘Evaluation of performance factors of fuel- 
oxidant mixtures,” Industrial and Engineering Chemistry, v. 43, Nov. 
1951, p. 2471-2475. 


This paper concerns itself with the solution of equations of condition 
for thermodynamic equilibrium in multicomponent gas mixtures. The re- 
sulting non-linear simultaneous equations are solved with the help of the 
Newton-Raphson method. The calculations were made on the IBM Card- 
Programmed Calculator,’ using the general purpose board on which the 
author collaborated. 

M. ABRAMOWITZ 
NBSCL 


1Sruart R. BRINKLEY, JR., G. L. WaGNER, & R. W. Situ, Jr., ‘General purpose 
ten-digit arithmetic on _ 1BM Cict Programmed Electronic Calculator,” Proceedings of 
Industrial Comput . New York, IBM, 1950. 


4. H. W. Brauau, W. G. SCHLINGER & B. H. Sace, “Evaluation of equa- 
tion of state constants with digital computers,” Industrial and Engi- 
neering Chemistry, v. 43, Nov. 1951, p. 2442-2446. 


In fitting the Benedict equation of state for the lighter hydrocarbons 
to experimental data, a technique of least squares is employed to evaluate 
the constants. The equations involved are complicated by a series of expo- 
nential terms and are arranged for solution by automatic digital computing 
machines. Tabular results for the Benedict equations for propane and 
methane, obtained with an IBM Type 604 Calculator, are presented. 


J. WEGSTEIN 
NBSCL 
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5. J. W. Donnett & R. W. Draper, “Absorption calculations by punch 
card calculators,” Industrial and Engineering Chemistry, v. 43, Nov. 
1951, p. 2449-2453. 


A procedure is given to eliminate the trial and error methods in accurate 
calculations made to determine the number of absorber trays and amount 
of absorber medium necessary to yield a given absorption. The systems of 
simultaneous equations obtained from the conditions of material balance 
and heat balance are solved by an iteration technique. Starting with pre- 
liminary estimates the iteration is carried through a number of cycles until 
convergence is obtained. The calculations are made on the IBM 602 Calcu- 
lating Punch and are so routinized that no judgment is necessary in the 
successive steps of the procedure. 


M. ABRAMOWITZ 
NBSCL 


6. R. M. Goopman, “A digital computer timing unit,” I.R.E. Proc., v. 39, 
Sept. 1951, p. 1051-1054. 


The EDVAC, a synchronous serial computer, requires a set of 48 timing 
pulses, each of which occurs at a particular pulse time within the 48-pulse 
minor cycle. These come from a: six-by-eight array of gates. Rows are 
governed by the tap points on a six pulse-time delay line which accepts 
every sixth clock pulse, and columns are governed by eight flip-flops which 
are triggered by the last timing pulse from each column. 


R. D. ELBourN 
NBS Electronic Computers 


7. GILBERT W. Kuro, ‘‘Monte Carlo method for solving diffusion problems,” 
Industrial and Engineering Chemistry, v. 43, Nov. 1951, p. 2475-2478. 


This article is a brief survey of some problems that have been handled 
by the Monte Carlo technique. It is pointed out that it has been customary 
but not necessary in solving problems in mathematical physics to formulate 
them as differential equations to be solved. Using Monte Carlo methods, 
highly realistic models representing very complex problems in classical 
physics or chemistry may be handled without going through the abstraction 
of the differential equation. The power and simplicity of the Monte Carlo 
method, when large scale computing facilities are available, is emphasized, 
and illustrations are given for several types of diffusion problems. 

J. H. Levin 
NBSCL 


8. FREDERICK J. MARTIN & Morris YACHTER, ‘Calculation of equilibrium 
gas compositions,’ Industrial and Engineering Chemistry, v. 43, Nov. 
1951, p. 2446-2449. 


The authors consider a set of non-linear, algebraic equations: 


DX = fi(xi, X2, Xn), 4=1,2,---,m 
j=1 
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which we will write 
(1) Ax = f(x). 


If x is the solution, and Z is an initial guess to the solution, the authors 
then wish to solve for the error 


(2) Ax 
by means of the following expansion: 
A(& + Ax) = f(@ + Ax 


or 
(3) [A — J(%)] Ax = f(a) — Az 


where J(Z) is the Jacobian matrix evaluated at Z. 

Equation (3) is a linear equation for Ax which is easily solved and the 
solution Ax is then substituted in (2) to yield the next approximation to x. 
The process is then repeated until the iteration converges. 

As a labor saving device, the authors suggest computing J(Z) just once 
and using the same matrix A — J(#) in (3) for all further iterations. 


LEON NEMEREVER 
NBSCL 


9. ASCHER OPLER & G. HeE!Tz, “‘Punched-card calculation of six- 
component distillation,’ Industrial and Engineering Chemistry, v. 43, 
Nov. 1951, p. 2465-2471. 


A procedure is given for the calculation of the plate composition of 
six-component continuous distillation in thermal balance. The equations 
resulting from the conditions of material balance and heat balance are 
solved by trial and error. The computations were made with the IBM 602A 
Calculating Punch. 


M. ABRAMOWITZ 
NBSCL 


10. ARTHUR Rose, R. Curtis JoHNson, & THEODORE J. WILLIAMs, ‘‘Step- 
wise plate-to-plate computation of batch distillation curves,”’ Industrial 
and Engineering Chemistry, v. 43, Nov. 1951, p. 2459-2464. 


The application of the IBM Card-Programmed Calculator to the pre- , 
diction of distillation operations is discussed. The problem is discussed for 
the case of non-ideal conditions such as carrying relative volatility, non- 
adiabaticity, varying plate holdup, beating effects, and condenser holdup. 
Comparison of the computed results with experiment is also made. 


M. ABRAMOWITZ 


11. ARTHUR Rose, R. J. LomBarpo, & THEODORE J. WILLIAMS, ‘Selective 
adsorption computations with digital computers,’’ Industrial and Engi- 
neering Chemistry, v. 43, Nov. 1951, p. 2454-2459. 


The authors discuss the approximate solution of the problem of adsorp- 
tion in liquids by considering the relations obtained from conditions of 
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successive material balance in different types of mixtures. The resulting 
system of equations is solved recursively on the IBM Card-Programmed 
Calculator using a general purpose board. 

M. ABRAMOWITZ 


12. ARTHUR RosE, THEODORE J. WiLLIaMs, & Harry A. Kaan, “Digital 
computers for trial and error calculations of distillation design,” Indus- 
trial and Engineering Chemistry, v. 43, Nov. 1951, p. 2502-2506. 


This paper discusses the application of the IBM Card-Programmed 
Calculator to distillation problems involving trial and error techniques for 
solution. Assuming certain basic conditions for the distillation, the authors 
determine the remaining product compositions employing the relations for 
material balance. A description of the machine procedure is given in some 
detail, and the results are described for the determination-of the reflex ratio 
under the assumption of constant relative volatility. The extensions to 
non-constant relative volatility and nonadiabatic operation are also dis- 
cussed. 

M. ABRAMOWITZ 


NEws 


Joint AIEE-IRE Computer Conference Committee.—This committee, with the par- 
ticipation of the Association for Computing Machinery, held a meeting in Philadelphia on 
December 10, 11, and 12, 1951. Some 900 persons attended the meeting. The program 
of the conference was the following: 


Monday, December 10 


Morning Session J. G. Brarnerp, Moore School of Engineer- 
ing, University of Pennsylvania, Chair- 


man 

Keynote Address W. H. MacWruiams, Bell Telephone Lab- 
oratories 

The UNIVAC System J. P. Eckert, Jr., J. R. Werner, H. F. 


WE su, and H. F. Jr., Eckert- 
Mauchly Computer Corporation, Division 
of Remington-Rand, Inc. 

Performance of the Census UNIVAC Sys- J. L. McPuHerson, Bureau of the Census, 


tem and S. N. ALEXANDER, National Bureau 
of Standards 
Afternoon Session I. Travis, Burroughs Adding Machine Com- 
pany, Chairman 
The Burroughs Laboratory Computer G. R. HoserG, Burroughs Adding Machine 
Company 
Inspection Trips Eckert-Mauchly Division of Remington- 


Rand, Technitrol Engineering Company, 
Burroughs Research Division, and Moore 
School of Electrical Engineering 
Evening Session H. E. Tompxins, Burroughs Adding Ma- 
chine Company, Chairman 
The Significance of Electronic Computers S. N. ALEXANDER, National Bureau of 
to Science and Management Standards 
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Tuesday, December 11 


Morning Session C. V. L. Smitrn, Office of Naval Research, 
Chairman 

IBM Card-Programmed Calculator J. W. SHELDON and L. Tatum, International 
Business Machines Corporation 

The ORDVAC R. E. MEAGHER and J. P. Nasu, University 
of Illinois 

Afternoon Session Dr. Mina REEs, Office of Naval Research, 
Chairman 

The ERA 1101 Computer F. C. MuLLaney, Engineering Research 
Associates 


The Operation of the Mark III Calculator G.E. Poorte, U.S. Naval Proving Grounds, 
Dahlgren, Virginia 
The University of Manchester Computing T. KiLBuRN, University of Manchester, 
Machine England, and B. W. PoLvarp, Ferranti, 
Ltd. 


Wednesday, December 12 


Morning Session J. M. Coomss, Engineering Research Asso- 
ciates, Chairman 

Whirlwind I Computer R. R. Everett and N. H. Taytor, Massa- 
chusetts Institute of Technology 

The EDSAC Computer M. V. WILKEs, University of Cambridge, 
England 

The National Bureau of Standards Eastern N. ALEXANDER and R. J. Stutz, National 

Automatic Computer (SEAC) Bureau of Standards 
Afternoon Session A. E. Samuz., International Business Ma- 


chines Corp., Chairman 
The Transistor as a Computer Component J. FELKER, Beli Telephone Laboratories 
Summary and Forecast J. W. Forrester, Massachusetts Institute 
of Technology 


Centre International de-Calcul Mécanique.—At a conference held in Paris in Novem- 
ber, 1951, the UNESCO voted to establish the Centre in Rome, Italy. This choice was 
influenced, no doubt, by the fact that Rome is the seat of the Istituto Nazionale per le 
Applicazioni del Calcolo, headed by M. Picone. A quarter of a century ago, Professor 
Picone established the Istituto di Calcolo in Naples. In 1932 it was moved to Rome, under 
its present name, and is located in the same building as the Italian National Council of 
Research. This building will be considerably enlarged to accommodate the Centre. It will 
contain at least one high-speed electronic calculating machine and will be staffed by spe- 
cialists recruited from all participating nations. The Centre’s yearly expenditure is expected 
to be around $130,000, which will be met by all the participating nations (about 20 in 
number) as well as by the Centre’s earnings from services rendered to private organizations 
in any country requesting them. On the other hand, these services will be available without 
charge to government scientists of the nations. 

A number of Italian scholars are already in the United States engaged in an extensive 
survey of electronic computation techniques and the engineering principles involved in the 
construction of high-speed calculators. Drs. Divo DaINELLI and Enzo APARO, members 
of Professor Picone’s Istituto, have been at the NBSCL for the past three months and have 
computed several important mathematical functions on the SEAC. They are planning 
an extended stay at the Harvard Computation Laboratory also, where they will join 
Gut1io RoptNo. The latter, an electronic engineer, replaced his colleague, MICHELE CANEPA, 
who returned to Italy after a year’s study of Harvard’s high-speed computers. 

The first post as American observer at the Centre is held by Dr. HERMAN GOLDSTINE 
of the Institute for Advanced Study, Princeton University. 
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UNIVAC Acceptance Tests.—On February 4-5, 1952, the second UNIVAC, con- 
structed by the Eckert-Mauchly Division of Remington Rand under NBS contract for the 
Office of the Air Comptroller, USAF, passed a magnetic tape reading and writing test which 
was the final test for its acceptance. The machine is now being moved from the factory in 
Philadelphia to the Pentagon Building in Washington, D. C., where its primary activity 
will be computing logistic programs. 

In this test UNIVAC read over 142 million decimal digits from magnetic tapes and 
wrote over 85 million on tapes in eight hours net running time and two hours down time. 
It ran without error or stoppage through 23 out of 32 fifteen-minute test units. In the other 
9 test units automatic checking circuits stopped the machine for 15 tape reading errors, 
for two malfunctions of tape driving mechanisms, for two tape defects, and for one tube 
failure in the computer. No errors escaped the automatic checking circuits. 

The same general test of computational ability that was given the first UNIVAC 
was also given to this machine. (See MTAC, v. 5, p. 176-7.) During the 19 twenty-minute 
test units, the only malfunctions were three tape reading errors which were detected by 
automatic checking circuits. 

The test of the Uniprinter required it to read from magnetic tape and to type 144,000 
characters among which occurred every typewriter symbol. This took four hours and 
fifty-four minutes running time plus nineteen minutes down time. Checking circuits stopped 
the printer six times, apparently because it picked up spurious pulses in the spaces between 
blocks of data on the tape. No errors were found in the copy. 

The Unityper test required 54,000 characters of instruction codes to be typed onto 
magnetic tapes and these to be readable by UNIVAC and by Uniprinter. Net typing time 
was 3} hours. UNIVAC read all the tapes correctly, but Unityper omitted characters at 
three points because tape slippage on Unityper had caused them to be recorded too closely 
together. 

In the Bureau of the Census’s year of experience with the first UNIVAC the computer 
has been remarkably reliable except for rather frequent tape reading error; therefore the 
tape reading and writing test for the second UNIVAC was revised to require nearly five 
times as much reading. For this reason the similarity in performance of the first and second 
UNIVACS in their tests really indicate noteworthy improvement in this respect. 
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13. Anon., “Aircraft and flight-control-system analog,” Instruments, v. 23, 
1950, p. 568, 570. 


A description with three pictures of the M. I. T. Flight Simulator con- 
structed under the supervision of A. C. HALL. This device consists of a 
continuous computer and a three gimbal flight table positioned by hydraulic 
servos. 

F. J. M. 


14. Anon., “Analog computer,” Instruments, v. 24, 1951, p. 772, one photo- 
graph. 


This article describes a differential analyzer suitable for solving systems 
of differential equations with constant coefficients constructed for M. I. T. 
by a Boston firm. The device is electrical, set up by means of patch cords, 
and can handle six equations in six unknowns. It “is not a commercial 
product.” 


F. J. M. 


| 
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15. ANnon., “Knows all angles,” Instruments, v. 24, 1951, p. 1178, one 
photograph. 

A drawing board instrument is illustrated and described which “gives 
true horizontal directions from oblique aerial photographs.” The device has 
three parts and is to be used with a specially designed slide rule calculator 
which gives certain required adjustments. 

F. J. M. 


16. Anon., ““USAF’s new jet instrument trainer,” Instruments, v. 23, 1950, 
p. 678, 680. 


Illustrated article on blind flying trainer (C 11 Link) containing a fixed 
continuous computer for the aerodynamic equations of flight. 


17. M. J. Bono, “Application of computing mechanisms to industrial 
instruments,” Instruments, v. 23, 1950, p. 614-616. 


The computing mechanisms are adding and multiplying linkages. By 
means of these one may construct devices to present directly a quantity 
computed from a number of measured quantities. 

F. J. M. 


18. B. CHANCE, F. C. WiLtiaMs, C. C. YANG, J. Busser, & J. HiGGins, 
“‘A quarter-square multiplier using a segmented parabolic characteristic,” 
Rev. Sci. Instruments, v. 22, 1951, p. 683-688. 


The multiplier described was developed for a fast analog computer. The 
parabolic characteristic is obtained by using 15 suitably biased diodes. 
The input to the squaring device must be positive. The squares of the sum 
and difference are obtained from the same squaring circuit on a time sharing 
principle, the amplitude of the square wave output is then proportional to 
the difference of these squares and is the required product. The parabolic 
characteristic covers a range of 0 to 25 volts within .4 percent, the output 
square wave has a frequency of 50 kilocycles per second and the overall 
accuracy was better than 1 percent. 

F. J. M. 


19. T. J. Conno.tiy, S. P. FRANKEL, & B. H. Sace, “Predicting phase 
behaviour with digital computers,” Elect. Eng., v. 70, 1951, p. 47. 


The authors discuss briefly the use of the I.B.M. 604 in certain thermo- 
dynamical calculations and refer to a more extensive monograph on this 
procedure by them, published by American Documentation Institute, 
Washington, D. C. 

F. J. M. 


20. J. T. CorTELyou, ‘Calculation of orifice meter charts by the square 
root planimeter method,” Instruments, v. 23, 1950, p. 1081. 


The planimeter mentioned but not completely described seems to be an 
ordinary planimeter suitable for use on the chart output of a recording 
device and the “square root”’ appellation refers to certain approximations 
made in using it. 


F. J. M. 
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21. M. J. Everitt, “A vernier method of using a slide rule,” Jn. Sci. 
Instruments, v. 27, 1950, p. 316-317. 


22. H. H. Goong, “‘Simulation—its place in system design,” I1.R.E., Proc., 
v. 39, 1951, p. 1501-1506. 


The author defines the simulation of a system as the “cut and try 
examination of its mathematical representation by means of a large analog 
or digital computer.’’ The design of a system is presented as progressing 
through a series of steps beginning with an analysis of the expected inputs 
and ‘a first order block diagram.’’ The author advocates that the early 
stages proceed by the use of linearized analysis but as the system takes on 
definite form he advocates that simulation in the above sense be used before 
large scale experimentation with actual systems. The latter is very expen- 
sive and should be limited as much as possible by “simulation.” For 
simulation purposes the author discusses among others the following topics: 
(1) Choice of Computer, (2) Programming, (3) Choice of Cases to be 
Treated, (4) Coding, (5) Machine Set-up, (6) Running, and (7) Data 
Reduction. In each case he presents a brief summary of experience and 
practice in the field. For instance in the choice of the computer, he lists 
types of systems for which automatic digital computers are most suitable 
and those for which continuous computers are better. However, the joint 
use of the two types of equipment is not emphasized. This paper is a valu- 
able summary of practice and experience in the field of simulation of systems. 

F. J. M. 


23. R. A. Harrincton, “Generation of functions by windup mechanisms,” 
Rev. Sci. Instruments, v. 22, 1951, p. 701-702. 


This instrument is based on the use of a suitably shaped cam on which 
a tape is partly wound. The tape also passes over pulleys and carries a pen. 
The displacement of the pen is then a prescribed function of the angle of 
the rotation of the cam. A method is described of cutting a cam for a given 
function and the limitations on the function are given. 
F. J. M. 


24. J. C. JAEGER, “A Schmidt mechanism for the approximate solution of 
the equation of linear flow in a medium whose thermal properties depend 
on the temperature,” Jn. Sci. Instruments, v. 27, 1950, p. 226-227. 


A drawing board instrument is described for the solution of the equation 


Ov ov 
(1) 2 (x2) = 


in which K, p, and ¢ are prescribed functions of v. Following van Dusen, 


the substitution @ = Ky f K dv, Ky = K(0) is made, so that (1) becomes 
0 


— = «(6)-'= where «(6) = K/pc. As in the derivation of the 


Schmidt method for the heat equation with constant properties, it is seen 
that if (2) is replaced by its usual finite difference approximation with 
(Ax)* = 2(At)«(@) then O(x,¢ + At) = $[0(x — Ax, t) + O(x + Ax, é)]. 
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In contrast with the customary Schmidt method for constant thermal 
properties, the required Ax depends here on 6(x,#). The instrument de- 
scribed produces the required Ax automatically so that as two points of the 
instrument are caused to trace the curve of 6(x, ¢) for a given t, a third point 
will draw the curve of 0(x, ¢ + At). The instrument utilizes a simple panto- 
graph linkage. For each function «(@) and choice of At, it is necessary to 
cut a plate cam representing (2(At)x(@))#. The various boundary conditions 
are treated as in the usual Schmidt process. 

A modification of the instrument is described for the solution of the 
problem in which a term A(v) is added to the left hand side of (1), repre- 
senting heat generation in the heat conduction problem. This requires the 
use of a second plate cam. 

J. H. WEINER 


Columbia University 
New York 27, New York 


25. H. JEFFERSON & L. M. Ericsson, ‘‘Mu-beta effect calculator,” I.R.E., 
Proc., v. 39, 1951, p. 1571. 


An additional method of using the Felker calculator. [I.R.E., Proc., v. 37, 
1949, p. 1204-1206; MTAC, v. 4, 1950, p. 175.] 


26. W. J. KENNEDY, “Chart reading by the Emco McGaughy integrator 
method,” Instruments, v. 23, 1950, p. 1082-1083. 


This article describes a mechanical computer involving square root cams 
and a disk integrator which utilizes instrument charts as inputs. 


F. J. M. 


27. E. Lakatos, ‘Checking analogue computer solutions,” I.R.E., Proc., 
v. 39, 1951, p. 1571. 


This letter indicates a method of checking an analogue computer against 
errors in scale factors and machine set up. For differential equations the 
machine values of the variables and their derivatives are substituted in the 
equations. 

F. J. M. 


28. F. J. LLEWELLYN, “A mechanical electrical unit for calculating structure 
amplitudes,” Jn. Sci. Instruments, v. 28, 1951, p. 229-230. 


This is a device to aid in the evaluation of multiple Fourier series. 
Multiplication of sin kx by sin ly is obtained by potentiometers. The re- 
solver action is obtained by a cycloidal arrangement of gears. One such 
unit is used to successively calculate the various terms. 

F. J. M. 


29. D. L. MacApaM, ‘Method of colorimetric integration using punched- 
card accounting machines,’ Optical Soc. Amer., Jn., v. 40, 1950, p. 
138-140. 
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30. J. H. Rospertson, “A simple machine capable of Fourier synthesis 
calculation,” Jn. Sci. Instruments, v. 27, 1950, p. 276-278. 


This machine is a Fourier synthesizer with a novel method of addition 
for the terms. The terms cause mirrors to rotate and, consequently, a beam 
of light successively reflected by all of them is displaced an amount approxi- 
mately proportional to the sum of the rotations. 

F. J. M. 


31. M. L. Unacar, “Probability paper,”” Permanent Records of Research 
and Development, Ministry of Supply, Shell Mex House, Strand, 
London WC 2, England (Sept. 1950). 


This monograph describes with samples three types of probability graph 
paper used in the British Ordnance Factories. Type (a) is based on the 
normal distribution function and can be used to test whether a distribution 
is normal and if it is, the mean and standard deviation. Type (b) is used 
for skew frequency distributions and Type (c) for a Poisson distribution. 

F. J. M. 


32. V. VaNnp, “A mechanical x-ray structure-factor calculating machine,” 
Jn. Sci. Instruments, v. 27, 1950, p. 257-261. 


The machine described is a mechanical device for summing trigono- 
metric series of the type developed by Kelvin. A method is described for 
varying the basic frequencies along a gradient in order to obtain the best 
fit to a given empirical curve. A new machine in which the summation 
process is by moments rather than by the Kelvin pulley method is being 
developed. 

F. J. M. 


33. J. F. Waters, “Instrument for measuring the slope of graphs,” Jn. 
Sci. Instruments, v. 28, 1951, p. 116. 


34. E. H. WINKLER, “Principle and design of a new type Stieltjes inte- 
grator,” Rev. Sci. Instruments, v. 22, 1951, p. 406-410. 
The objective of this device is to evaluate 


f F(x) dH (x). 


F(x) and H(x) are given as curves and a mechanical device plots the curve, 
given parametrically by v = F(x), « = H(x). The area under this curve 
corresponding to the given integral is then obtained by a planimeter. An 
electric follower is described which traces a curve drawn with conducting 
ink. The tracer is to remain on one side of the curve. 

F. J. M. 


NOTES 


132. FURTHER STATISTICS ON THE DiGits oF e.—An account of the 
calculation of the first 2500 digits of e on the ENIAC appears in MTAC, 
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v. 4, p. 11-15 and p. 109-111. In December, 1950 GrorGE REITWIESNER 
of the Ballistic Research Laboratories at Aberdeen furnished the Computing 
Service with a deck of punched cards bearing the remainders from the 
ENIAC calculation at 2520, 2538, 2556, 2574 and 2592 D. The remainders 
at 2520, 2538 and 2574 D were punched on new cards and the calculation 
was extended on the 602A. Originally the machine was wired to calculate 
5D of quotient and 4D of new remainders. Later the process was speeded 
up to punch 12D of quotient and 4D of remainder. The job ran as a fill-in 
project for the 602A until the three runs had all passed 3000D. A new set 
of remainders, at 2556D, was then punched and the calculations were re- 
peated on the 604. The wiring, to calculate 20 digits of quotient in one pass 
through the 604, as well as the entire calculation to check, was done by 
ORVILLE MARLOWE. Thus all digits from 2557 to 3002 were done four 
times. The remainders at 2945, 3005, 3010 and 3023D are available should 
anyone wish to extend the calculation. 

The x? analysis of the digits in each hundred (MTAC, v. 4, p. 110) may 
be extended as follows: 


n A, = Xa* Pn 
2600 2.04 .0092 
2700 2.54 .0202 
2800 2.53 .0199 
2900 2.67 .0240 
3000 2.95 .0339 


FRED GRUENBERGER 
Computing Service 
University of Wisconsin 


133. On FINDING THE ELECTRONIC CYCLE CAPACITY OF THE IBM Type 
604 CompuTER.—Frequently, when programming problems on the IBM 
Type 604 Electronic Calculating Punch, it is convenient to know the capacity 
of the machine in electronic cycles. This is especially true when one is 
attempting to partition a long empirical formula so as to effect the solution 
with a minimum number of runs on the machine. The user of the machine 
does not ordinarily have access to the driving oscillator so as to be able to 
change the operating frequency. Furthermore, extended operation at other 
than the prescribed frequency may be detrimental to the machine. Hence, 
the scheme described herein was devised for finding the capacity of the 
machine at the frequency which has been set internally. According to the 
operating manual (p. 65), the capacity is 252 electronic cycles. This, how- 
ever, will vary slightly from one machine to another and within a given 
machine over a period of time. The following is a method for determining 
the machine capacity in 5 to 10 minutes. 

We start by key punching a card with all 9’s in columns 1-12 and a 7 
in column 13. The 521 is wired to read columns 1-8 into FS 3-4 and columns 
9-13 into MQ, to emit O into FS2, and to stop in the case of an unfinished 
program. In each of the first four program steps ROFS3-4, MULT +. This 
will require (48)(4) = 192 cycles. On program step 5 ROFS2, RIMQ (1 
cycle) and finally on program step 6 ROFS3-4, MULT + (5 cycles). The re- 
maining programs are unwired, hence in the case of a 60 program step 604 we 


= 3 


5 
ot I 
7 
ie 

E 
F 

] 


QUERIES 125 


have 54 additional cycles or a grand total of 252 cycles. By varying the 
number emitted into FS2, the number of cycles in program step 6, and 
hence in the total, may be increased. By varying the number read into MQ 
from the card the number of cycles may be increased or decreased by 
multiples of 4. If the machine has fewer cycles the unfinished program will 
cause it to stop. Usually key punching 3 or 4 cards and moving one wire 
on the 521 control panel 2 or 3 times will determine the capacity of the 
machine to the nearest cycle. 

If the user has a 20 or 40 program step 604 the programming must call 
for 40 or 20 electronic cycles less respectively. This can be accomplished by 
simply changing the number read into MQ. 


Anpress O. RIpGway 
5618 Chillum Heights Drive 
Hyattsville, Maryland 


134. Fritz Empe.—Fritz Emde, professor emeritus of the Technische 
Hochschule, Stuttgart, died there on June 30, 1951; he would have been 
78 years old on July 13. In MTAC, v. 3, we published (p. 397-398) some 
biographical notes, and (opp. p. 333) a portrait. 


R. C. ARCHIBALD 
Brown University 


Providence, R. I. 


135. RALPH ERNEST PowErS.—This amateur mathematician died on Jan. 
31, 1952, at Puente, California. He would have been 77 years old on April 27. 
Mr. Powers was more responsible than any other man for the demonstration 
of the failure of Mersenne’s conjecture. He proved that 2° — 1 and 2 — 1 
were primes, and that several other Mersenne numbers were composite by 
long and laborious desk machine calculations. He was not aware of the 
discovery, the night before his death, of two new Mersenne primes (MTAC, 
v. 6, p. 61). Mr. Powers was born in Fountain, Colorado, and spent most 
of his life in Denver. 

D. H. L. 


QUERIES 
41. A DEFINITE INTEGRAL.—Does the integral 


f exp (—su — u-~*) du 
0 


have a generalized asymptotic form for large positive values of z, say 
z > 10? The integral is a solution of the differential equation 


+ + 2y = 0, 
Leo A. AROIAN 
Hughes Aircraft Co. 
Culver City, Calif. 
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CoORRIGENDA 


CORRIGENDA 


. 9,1. — 1, for up read only up. 

. 16,1. — 18, for ka read $kr. 

. 16, 1. — 8, for 2kw read kr. 

. 17, footnote 9, for 155 read 255. 

. 26, 1. 15, for + read + y. 

. 31, 1. 5, for y read y. 

. 32,1. — 9, for 956 read 957. 

. 41, 1. 8, insert authors G. FLAKE, A. JONES. 
. 41,1. 15, for = read p =. 


V. 6, p 
V. 6, p 
V. 6, p 
V. 6, p 
V. 6, p 
V. 6, p 
V. 6, p 
V. 6, p 
V. 6, p 
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